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ABSTRACT 


Thermal design of solar receivers is commonly accomplished via 
approximate models, where the receiver is treated as an isothermal box 
with lumped quantities of heat losses to the surroundings by radiation, 
conduction and convection. 

These approximate models, though adequate for preliminary design 
purposes, are not detailed enough to distinguish between different 
receiver designs, or to predict transient performance under variable 
solar flux, ambient temperatures, etc. A computer code has been written 
for this purpose and is given the name "HEAP", an acronym for Heat Energy 
Analysis Program. HEAP has a basic structure that fits a general heat 
transfer problem, but with specific features that are custom-made for 
solar receivers. The code is written in MBASIC computer language. 

This document explains the detailed methodology followed in 
solving the heat transfer problem, and includes a program flow chart, an 
explanation of input and output tables, and an example of the simulation 
of a cavity-type solar receiver. 
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SECTION 1 


INTRODUCTION 


1.1 BACKGROUND 

The need for an accurate prediction of the thermal behavior 
of solar receivers, to be used as a basic tool for their performance 
comparison, has become increasingly important with the rapid develop- 
ment of solar energy technology. A sufficiently general analytical 
model is needed which is both an adequate idealization of the physical 
system and is capable of expressing itself by simple mathematical descrip- 
tion. 


Thermal design of solar receivers is commonly accomplished 
with approximate models where the receiver is treated as an isothermal 
box with lumped quantities of heat losses by radiation, conduction and 
convection to the surroundings. These approximate models, though adequ- 
ate for preliminary design purposes, are not detailed enough to distin- 
guish between different receiver designs, to predict transient perfor- 
mance, or to accurately express the performance sensitivity to variations 
occurring in the system parameters. 

Since the solar receiver is part of the energy collection 
subsystem, and located in a very strategic spot in between the con- 
centrator and the energy conversion subsystem, its thermal performance 
and concept selection affect the system installation and operation cost. 

A detailed model is ; therefore, needed not only to predict the per- 
formance under varying solar flux, ambient temperatures, and local heat 
transfer rates, but also to detect locations of hot spots and metalluri- 
cal difficulties, and to predict the performance sensitivity to neighbor- 
ing component parameters. 
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Hie computer code written for this purpose is presented in 
this report and is given the name HEAP, an acronym for Heat Energy 
Analysis Program. HEAP has a basic structure that fits any general heat 
transfer problem, but with specific features that are "custom-made" for 
solar receivers. The processes that led to the writing of this new 
code and the criteria to be satisfied are explained in the following 
sections. 


1.2 SURVEY OF OTHER AVAILABLE PROGRAMS 

Before constructing the HEAP code, a computer library search 
was conducted to review current programs and select one that could be 
applied. Five computer programs were identified in the search: 

(1) LOHARP: Lockheed Orbital Heat Rate Package, contains 

three subprograms for calculation of radiation view 
factors, radiation constant, and radiation orbital 
heat rate. The program is applicable to outer space 
related problems. The program is documented in 
Reference 1. 

(2) TAS: Thermal Analysis System, a JPL-developed pro- 

gram for steady state problems only, written in 
FORTRAN IV. No convection heat transfer or fluid flow 
across nodee are considered (see Reference 2) . 

(3) THERM: Therma l Analyzer Computer Code is another code 

in FORTRAN IV, which is capable of handling 1000 nodes 
of a problem with steady state or transient heat 
transfer. Available documentation (Reference 3) is 
not complete or clear in its assumptions and limita- 
tions . 

(4) SINDA: Systems Improred Numerical Differencing 

Analyzer is a code expanded by TRW from a version 
written by Chrysler Corporation. It is a general, 
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versatile but expensive code which handles both 
steady and transient problems. Users may add 
their own subroutines after becoming familiar with 
its massive documentation (Reference 4) . 

(5) HEAT: HEAT Transfer Computing Program, developed 

by W. A. Beckman at the University of Wisconsin. 

HEAT is basically the same as TAS except that 
transient solutions are added. HEAT is written in 
FORTRAN V language (see Reference 5) . 

1.3 SELECTION CRITERIA 

Before selecting the computer ptogram most suitable for 
advanced receiver studies, some criteria were set to bound the selection 
process of any one of the above computer codes: 

(1) The computer program should include analyses of all 
modes of heat transfer by conduction, convection, 
longwave (infrared) radiation, solar radiation and 
the convective energy carried in and out the receiver 
by flowing fluids. The program should also allow 
for possible heat generation or loss due to chemical 
reactions, electrical resistances, etc. 

(2) The computer program must handle both steady state 
(or quasi-steady state) and transient problems. This 
requirement is essential in testing receivers of 
different response times under variable daily solar 
flux, ambient temperatures, or flow rates. 

(3) For ease of use, the computer program should be a 
self-sufficient and complete package with minimum or 
no dependence on other computer program outputs. This 
means that the user should not be required to put 
together the input data for one program using sub- 
routines of other programs as a prerequisite or to 
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reach the final answers of the problem by working 
in more than one program simultaneously. 

(4) The program documentation should be dear, staple, 

and self explanatory to shorten the user's start-up time 
and should explicitly show its flexibility and limita- 
tions . 

(5) The computer program should be closely associated with 
this particular application to solar receivers, but 
generally expressed to handle many possible receiver 
configurations. This criterion saves time consumed 

in running complicated, general programs that have too 
many variables irrelevant to the problem under study. 

The five computer codes previously identified were evaluated 
against these selection criteria. The first program, LOHARP, was 
dropped because of the fifth criterion since it handles mostly outer 
space objects with no convection heat transfer. The second program, 

TAS, was also dropped because of the first and second criteria since 
it handles steady state problems only, with no heat or flow convection. 

The third program, THERM, was omitted for lack of good documentation, 
required by the fourth criterion. The fourth and fifth programs, SINDA 
and HEAT, though comprehensive at times, both lack the completion speci- 
fied by the third criterion. They both require radiation view factor 
data which the user may either calculate or obtain from the output of 
the LOHARP program, for example. HEAT assumes that the user will pro- 
vide all thermal and fluid conductances, which can be quite a tedious 
process. For instance, for a 30-node problem, at least 450 conductances 
are needed. The SINDA program was further downgraded for being too 
complex and still needs extra subroutines to be added. The fifth, and 
the only hopeful program, HEAT, was then carefully investigated to see 
how additional receiver subroutines could be tied to it. Unfortunately, 
the HEAT structure and assumptions were found Inadequate in handling 
many computational procedures. For the Euler's minimum time step of 
the trans'jnt solution, HEAT approximates the radiation effect, which 
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is dominant in solar receivers, as equal to one-half of all other effects, 
an approximation which may not be adequate. Also, since the user pro- 
vides all conductances, HEAT makes no distinction between different 
surface areas and neighboring nodes. For each node, only one surface 
area for all modes of heat transfer is used in HEAT, which is not adequ- 
ate for a detailed heat transfer analysis. 

Finally, a decision was made to save time and effort con- 
sumed in correcting, adding, or deleting in the candidate programs to 
suit solar receivers, by writing a new computer code that satisfies the 
above criteria and fits the sp ^ial needs. The following sections of 
this report explain the detailed methodology followed in solving the 
heat transfer problem and show the flow chart of the program. A solved 
example is also provided in Section 4 to illustrate the program use for 
a cavity-type solar receiver. 
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SECTION 2 


PROGRAM METHODOLOGY 


2.1 GENERAL LAYOUT 

The methodology used In the heat transfer calculations is a 
well-known numerical technique which has enjoyed wide application since 
the advent of high-speed digital computers. The system (receiver) under 
consideration is discretized in space by nodes as shown in Figure 2-1. 
Heat exchange between any two arbitrary nodes (i) and (j) is illustrated 
in Figure 2-2. The heat transfer rate equations are written with all 
modes of heat transfer included. The net energy stored in each node 
(sometimes called the energy "residual") is calculated using th« lirst 
law of thermodynamics, including the energy exchange to and from each 
node as a result of neighboring nodes. 

Spatial nodes are divided into two categories: a) equilibrium 
nodes (or finite heat reservoirs) whose temperature is dependent on the 
net energy stored; and b) source or sink nodes (or infinite heat- 
reservoirs) which gain or lose any amount of heat without changing 
temperature. The program structure differentiates between these node 
categories for both steady state and transient solutions. Another 
advantage was taken from the axisymmetric geometry of receiver nodes 
by adding a special radiation view factor subroutine. To avoid 
excessive precalculations of thermal conductances, the program is built 
to compute these conductances with minimum information requested from 
the user. 


For the steady state solution, the energy residual to each 
equilibrium node must be zero. A Newton-Raphson Iteration solution is 
used to solve for: a) the temperature distribution of equilibrium 

nodes; and b) the residuals of source/sink nodes. 
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Figure 2-2. Heat Exchange Between Nodes (i) and (j) 
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For the transient solution, the "forward" finite difference 
numerical technique is used. By this technique, the new temperature 
vector at time (t + At) is calculated using the old temperature vector 
at time (t) a.ic all other parameters evaluated at time (t). To ensure 
stable computations (convergence) a special formula is used for the 
critical time step (At) based on the methodology stated in Reference 7. 

To minimize computational time consumed due to nodes having a very small 
time step (At), the concept of zero-capacity nodes was included in the 
transient solution. 

2.2 BASIC ASSUMPTIONS 

The following assumptions and idealizations were made in the 
heat transfer calculations: 

(1) The physical system under study (solar receiver) is 
treated as a collection of a number of isothermal 
nodes with uniform optical and physical properties 
throughout each node. Nodes are preferred to be of 
axisyrametric geometry (plane discs, plane rings, 
cylindrical discs, cylindrical rings etc.) to make use 
of the special radiation view factor subroutine. The 
axis of symmetry must coincide with the receiver axis. 

The condition of axial symmetry is preferred, but is 
not mandatory. 

(2) The nodal specific heat, thermal conductivity, and density 
are assumed independent of operating temperature. The 
properties may be taken at an estimated average work- 
ing temperature as an approximation. Also, nodal 
infrared and solar emittances are assumed independent 

of both directivity and operating temperature. A 
second phase of the program would consider the above 
variations if they showed serious effects on the 
computational accuracy. 
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(3) The solar flux incident on each node is assumed 
uniform. Therefore, the user should be careful when 
dividing the receiver into nodes to avoid distorting 
the flux profile inside or outside the receiver body. 

(4) The radiation emitted by a node surface is assumed to 
be diffuse. In general, the radiation reflected from 
a node surface can be diffuse, specular or a combina- 
tion of diffuse and specular parts. In the present 
model, only diffuse reflection was assumed, with a 
zero specular component. 

(5) The optical properties of a node surface are assumed 
to be divided into two bands in the electromagnetic 
spectrum. This "semi gray" assumption allows separate 
calculations of radiation energy exchange in the short 
wave (solar) band and the longwave (infrared) band. 
Consequently, the radiation energy in the program is 
allowed to be absorbed, reflected or emitted in the 
longwave region, and is only absorbed or reflected 

in the short wave region. 

(6) The nodal internal energy generation or dissipation 
by electrical, chemical or nuclear effects is assumed 
uniform over the node surface and independent of the 
temperature. 

Based on the above assumptions and idealizations, the heat transfer rate 
ouations are written and analyzed as described in the following sections. 

2.3 CONDUCTION HEAT TRANSFER 

For two neighboring nodes (i) and (j), exchanging heat by 
conduction across the interface area A(i,j), the heat transfer Q( i ) con ^ 
is given by 
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[T( j) - T( i) ] 


( 1 ) 


where 


O^cond ' £ C(1 -J>cond 




cond 


A(i,J) 

Mil + MU 
K(i) K(j)J 


( 2 ) 


K is the thermal conductivity, and L is the conduction distance between 
the center of each node and the separating surface. 


The value of C(i,j) could be a function of time and/or 
temperature or any other variable. However, the conductance C(i,j) is 
assumed constant and uniform throughout the node. Because of the 
symmetric nature of the conductance matrix, only half the input data 
are needed for A(i,j), where the half-cut is made along the diagonal. 


2.4 CONVECTION HEAT TRANSFER 


The convection heat transfer across the boundary layer of 
solid-liquid, liquid-gas or solid-gas interfaces, can be put in the 
following form between two neighboring nodes (i) and (j) of different 
phase: 


«<«co„v ■ f C(i -J>conv • 1 T <3> - «») < 3 > 

where the conductance C(i,j) „ , is given by: 

conv 

CU.jJconv " h(i ’ j) ’ A(i,j) (4) 

and where h(i,j) is the convective heat transfer coefficient. 

C(i,j) is treated in a manner that allows its value to be changed 
,J conv 

with time. This is made simply because the convection coefficient h(i,j) 
is a dominant function of the fluid mass flow rate which may not be con- 
stant during the receiver operation. Furthermore, because of the 
symmetric nature of the h(i,j) matrix, only half of it (cut along the 
diagonal) is needed as input data. 


2-6 


2.5 


FLUID FLOW AND ENERGY EXCHANGE 


Consider a node (i) surrounded by two adjacent nodes (j) and 
(k) with a fluid flowing at a rate m f from (s) to (i) to (k) as shown 
in Figure 2-3. An energy balance on node (i) taking into account the 
energy entering at T(s) and leaving at T(i) would yield 

Q(i) f - m f Cp f (i) [T(s) - T(i) ] (5) 

which will be written as 

Q(i) f = C(i,s) f [T(s) - T(i) ] (6) 

or 

Q(i) f - £ C(i,j) [T( j) - T(i) ] (7) 

j 

where T(j) is the temperature of all neighboring nodes in direct con- 
tact with node (i) . For a single fluid stream crossing node (i), there 
will be only two surrounding nodes, one upstream and the other down- 
stream. The flow conductance C(i,j) will then be taken as [m^dJCp^d) ] 
if the flow goes from (j) to (i), and zero if the flow goes from (i) 
to (j). The matrix C(i,j) f will then be asymmetric such that C(i,j) j 4 
C( j , i) and all elements are needed as input data. For boundary nodes 
that present the inlet fluid zone or outlet fluid zone as shown in 
Figure 2-4 the next neighboring node in the closed loop should always 
be considered. In the program, one input number (+1) will be entered 
in the [MASFLO] matrix representing the direction of the flow to each 
node. Care must be taken when these boundary nodes are treated as 
source/sink nodes and when dealing with more than one fluid loop inside 
the receiver. 
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FLOW 

m f 



(a) Fluid Flow Acrou Node (I) 


TEMPERATURE 



I 



Figure 2-3. Fluid Convection Between Nodes 
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Figure 2-4. Closed Sequence of Fluid Nodes in Receiver Loop 
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2.6 


INFRARED (IR) AND SOLAR RADIATION 


The analysis of heat transfer by radiation Is divided Into 
two wavelength regions: a) an infrared or long wavelength region in 

which infrared properties apply; and b) solar (or short wavelength) 
region in which solar radiation properties apply (a wavelength region from 
0.3-5u). The solution methodology is: 

For an enclosure composed of N surfaces as shown ip 

th 

Figure 2-5, the radiosity (J) streaming out of the i node is given 
by 


J(i) - X(i) + p(i)G(i) (8) 

where p(i) is the gray reflectance of the i th node, G(i) is the 
irradiation or flux streaming onto the node and X(i) is the thermal 
excitation at the i t * 1 node. For infrared or longwave radiation, the 
excitation (X) can be written as 

X(i) » oe(i) T 4 (i) (9) 

where a is the £ ttf an-Boltzman radiation constant. 

On the other hand, the relationship between the radiosity 

[J] and irradiation [G] for each node can be determined using energy 
conservation principle as 

A(i)G(i) - A(l)F(l,i) J(l) + A(2)F(2,i) J(2) + 

A(K)F(K,i) J(K) + ... + A(N)F(N,i)J(N) (10) 

where A's represent surface areas and F(k,i) is the view factor 
(sometimes called geometric, configuration or form factor) from node 

(K) to node (i). Using the view factors reciprocity relation pre- 
sented by 


A(K) F (K, i) - A(i) F (i,K) 


(ID 
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NOOf 0 


RADIO SITY J (i) IRRADIATION G 0) 



NODE N 


NODE0 



NODE (!) 


Figure 2-5. Raliosity and Irradiation in Enclosure 

Equation (10) is then reduced to the form 

G(i) - FU.DJll) + F(i,2)J(2) + ... 

+ F(i,K) J(K) + ... + F(i,N) J(N) (12) 

or in the matrix form to 
[G] - [F] [J] 

Equation (12) can be also put in the indicial form 

G(i) - £F(i,k) • J(k) (13) 

k 

Combining Equations (8) and (13), we can write 

J(i) - X(i) + P(i) [F(i,l)J(l) + F(i,2)J(2) + ...) (14) 
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Using Equation (9), the net radiant flux into i C ^ node is expressed as 

Q"(i) - [1- p(i) ] G(i)-X(i) (21) 

Combining with Equation (18) written in the form 

G(i) =LF(i,K)X(K) 
k 

then 

Q"(i) - 

Equation (22) can be applied to both the solar band and IR band of 
radiation. These cases are discussed below. 


E(l- p(i) 3 F (i.k) . X(k) 


-X(i) 


( 22 ) 


2.0.1 Energy Exchange in the IR Band 


The excitation X(i) is given by Equation (9). Accordingly, 
the net radiant flux into the i C ^ node becomes 





e(i) • F(i,k) • e(k) • T 4 (k) 


- e(i) T 4 (i) 


(23) 


In a uniform temperature enclosure, the net radiant flux Q"(i) is zero. 
Therefore, 


T(l) - T(2) - T(3) - ... » T(N) - T 


For example, for the first node in this case 


- e(l) (F(l,l) e(l) + F( 1 , 2) e(2) + ...] - e(l) - 0 
oT 


or, in general terms 


e(i) E^U.k) e(k) - e(i) 
k 

Pd£CEDING PAGE BLANK NOT FILU&k 1^1 
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The left hand side is sometimes called Hottel's F-script matrix, J'd.k), 
where 

E ^F(i.k) - e(i) ) 

k j (24) 

Z F (i,k) e (k) - 1 

k ] 

Equation (24) can still hold for a non-uniform temperature enclosure and 
by combining Equations (23) and (24) , then the IR radiation energy exchange 
to node (i) of area A(i) is 

Q(i) TR - oA(i) • E e (i) . F(i,k) • e (k) • [T 4 (k) - T 4 (i)] (25) 

IR k 

2.6.2 Energy Exchange in the Solar Radiation Band 


For the solar band, the IR properties p, e, are replaced 
by p*, e* for differentiation. The solar excitation X* on node (i) in 
this case is as shown in Figure 2-6 and is given by 

X(i) = p(i) • Flux (i) (26) 

Equations (15), (17), (18), and (19) still apply, but use the optical 
properties with an asterisk for the solar spectrum. 


A*(i,k) 

= 6 (i,k) - p *(i)F(i,k) 

(27) 

[J*] * 

[A*]" 1 • IX*] 

(28) 

[G*] = 

[F*J . [X*] 

(29) 

IF*] - 

[F] • [A*]' 1 

(30) 
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•K0 -p(0 G(l) *pO) FUiX(l) 



Figure 2-6. Radlosity and Irradiation in the Solar Band 


The net solar radiation flux retaining in the i C ^ node is expressed as 

Q"* (i) - G* (i) - p *(i) G*(i) + Flux(i) - p*(i)Flux(i) (31) 
Q"* (i) " [1 - p*(i) jG*(i) + o*(i)Flux(i) 


Using Equation (29), then 

Q"*(i) e*(i)F*(i»k)X(k) + e*(i)Flux (i) 

k 


From Equation (26), the net solar energy exchange on node (1), Q (i) is 


QU) S 


c*(i) *F*(i,k) *p(k) •Flux(k) + e(i)Flux(i) 
k 


A(i) 


s 

(32) 


where 


F* (i.k) 


£F(i,s) 


*-l 

A ( 8 , k? 


(33) 
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From the above analysis , Equation (25) is the final form 
to be used for the longwave (infrared) radiation exchange, and Equation 
(32) is the final form to be used for shortwave (solar) radiation 
exchange. 


2.7 NODAL ENERGY BALANCE 

For a particular receiver node (i), the energy is exchanged 

with the neighboring nodes by many modes of heat transfer such £<s: 

1) longwave (infrared) radiation, Q( i) IR J 2) shortwave (solar) radiation, 

Q(i) ; 3) conduction, Q(i) 4) thermal convection, Q(i) ; 5) 

s cona ' conv 

fluid convection, Q( i) ^ ; and 6) internal generation (or dissipation) by 

electrical, chemical or nuclear effects, Q(i) . The final mathematical 

P 

form of each heat transfer mode is written as follows: 

The net longwave (infrared) radiation received by the i C ^ 
node from all neighboring nodes is given by 

Q(i) IR 'L c (i.j) IR [T 4 (j> - T 4 (i)] (34) 

j 

where T is the absolute temperature and C(i,j)^ is the IR radiation 
conductance between nodes (i) and (j) . can be derived from 

Equation ( 25 ) . 

The final forms of the shortwave (solar) radiation, Q(i) g 
and the internal energy generation Q(i)^ are assumed independent of 
operating temperature and are constants. The energy transferred by 
conduction to the i^ node from neighboring nodes is given by Equation 
(1) as 

‘K 1 >co»d-S C<1 'J'co„d |I<J) - 1(1,1 

Also, for thermal convection, the ret energy convected to 
the i node is given bv Equation (3) stated before as 

Q(i) »E(C(i.j> [T( j ) - T(i) ] 

conv conv 

where C(i,1) is the thermal conductance for convection. 

J conv 
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th 

For fluid convection, the heat extracted or dumped to the i n node 
is given by Equation (7) as 

Q(i) f - £C(i,j) f [T(j) - T(i) ] 
j 

where C(i,j>£ is the fluid conductance for flowing fluid. 

The first law of thermodynamics applied to the i^ node gives: 

DN(1) V(l) CPU) - Res(i) - 

^IR + Second + "‘“conv + «<«f + «« s *»<«p <’» 

where Res(i), DN(i), V(i) are the net energy stored, density and 

f h 

volume of the i node, respectively, and At is the elapsed time 
increment . 


2.8 MATRIX FORMATS 

To facilitate the manipulation processes, Equations (1), 
(3),, (7), and (34) are put in matrix or vector forms. Generally, 

Q(i) -£c(i,J) [T n (j) - T°( i) ] (36) 

j 


where n is an exponent that can be either 1 or 4, and C(i,j) is a 
general indiclal form for conductance matrix. Equation (36) can be 
written in the expanded form 

Q(i) - C(i,l) [T n (l)-T n (i)J + C(i,2) [T n (2) - T n (i)]+ .. 
or 

Q(i) - C(i,l) T n (l)+C(i,2)T n (2)+. .- [C(i,l) + C(i,2)+ C(i,3)+. . ]T n (i) 


2-17 


For the equilibrium nodes, the set of non-linear energy 
equations are solved by iteration starting from an initial estimate of 
the temperature vector. The residuals will be equal to zero at steady 
state. The most convenient method used is the iterative Newton— Raphson 
method. The method linearizes Equation (39) in the region of the initial 
temperature estimate using the Taylor series. 

Therefore, 

6Res(i) = E jf f jp 5 T(j) (41) 

or 

[6Res] « [JAC] • [<5T] 

6Res(i) = E JAC ( i , j ) • <$T(j) (42) 

j 

where [JACj is the Jacobian matrix, given by using Equations (40) and (41) 
as 


JAC(i,.D - = 4C(i,j) TB T ( j ) 3 + C(i>j) eon( j + C(itj) 


IR 


cond 


conv 


+ ( i > j ) 


(43) 


Equation (42) can be written as: 


(K+l; (k) (K) 

[Res(i) - Res(i) ] - JAC(i,j) 


(K+l) (K) 

T( j) -T(j) 


or 


(K+l) 

[Res] - [Res] 


(K) 


[JAC] 


[T ] (k+1) -[t] ( k) * 


(44) 
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where K is the number of iteration process. Equation 44 can also be 
written as 

(K+l) (K) -1 (K) 

[T] = [T] - [JAC] * [Res] (45) 

Reference 6 states that six iterations or less are usually 
sufficient when initial estimates are far from actual steady state 
values. An exceptionally good initial guess does not necessarily save 
computer time or cost since it may eliminate only one or two iterations. 
Therefore, the user should not put much effort into the initial tempera- 
ture estimate accuracy. 


2.10 TRANSIENT SOLUTION 

After the elapse of a small time interval Ax, Equation (35) 
gives the slope of equilibrium nodes temperature array T as an approxi- 
mation of the derivative 3T/3x using the finite difference. The 
"forward" finite difference method will be used as it is the most 
adequate one The new nodal temperature T (i) at (x + Ax) is calcu- 
lated using the old temperature T (i) and all other properties at the 
old time (x) , i.e. , 

(x+Ax) (x) (x) 

" I(i) + M(lfcRi) Res<1> < 46) 

where M(i), Cp(i) and Res(i) are the mass, specific heat and energy 
residue at the i c node. The time step (Ax) is subject to computation 
convergence and stability criteria. To use the simple method applied 
in Reference 7 the residual Res(i) is written in the following 
linearized form 

Res(i) - £D(i,j) [T(j) - T(i) ] + E(i) (47) 

j 
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where 


DU.j) - «!.«„ [ 4? I ^‘-J’cond + C(1,j) conv + Cii ’»t 


E(i) » Q g (i) + Q p (i)*0 


combining Equations (47) and (48) , then 


(48) 


(t + Ax) 
T(i) 


(t) 

T(i) 



At 

M(i)Cp(i) 



ED(i,j)j +^ED(i,j) T(j) + E(i)j 

(49) 


Since all the coefficients D(i,j) and E(i) are always positive or at 
least zero, the coefficient of T(i) should be selected to be either 
positive or at least zero in order not to violate the thermodynamic 
principles . Accordingly 


At(1) 

min 


£ M(i)Cp(ij 
D(i,j) 


j 


(50) 


the denominator £D(i,j) can be expressed in terms of the modified 
i 


matrices [C] in Equation (38) as: 


£ D(i,J) 
j 


£c(i,j) 

1 3 

+ C(i, j) 


IR 


conv 


T( j) - T(i) 
+ C(i,j) f 


+ 


£c.(i,j) cond 


Eo(i,j) 

j 




LC(i,i) IR + C(i,i) 


cond 


+ C(i,i) 


conv 


+ C(i, 




(51) 
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where [LC] is the linearized 1R radiation conductance given by 



The matrices [LC] and [LC] are related by the same expression, 

Equation (38). For each node, Equation (50) is used to find the minimum 
time increment At (i) for stable computations of the i t * 1 node. A 
scanning process is then made to find the minimum of these minima to be 
used for the iteration process for all nodes, excluding source/sink 
nodes . 


2.10.1 Zero-Capacity Nodes 

Some nodes have a very small time increment, At (i), which 
reduces the overall time increment At to such a small value that the 
computer execution time and cost increase dramatically. These types 
of nodes have either a very small thermal capacity [M(i).Cp(i)] or a 
small thermal resistance and they do not store a significant amount of 
energy during the exchange with the surrounding nodes. The transient 
analysis in the program considers a minimum time increment (MINDTU) 
below which such nodes can be treated as zero-capacity nodes. 

Zero capacity nodes are momentarily in thermal equilibrium 
with neighboring nodes with no energy storage during intervals of At. 
The temperature vector at time t + At Is determined by the Newton- 
Raphson iteration process similar to the one described under steady 
state (Section 2.9) using the old neighboring node temperature at time 
t. For other equilibrium nodes that are not ’’zero-capacity” nodes, the 
new temperature vector at time t + At is determined by Equation (46) 
baaed on the old temperature vector at time t. In both cases, source/ 
sink nodes are not affecting the minimum time increment selection. 

The detailed flow chart of calculation is presented in Appendix A. 
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2.11 


INTERPOLATION SUBROUTINE 


This subroutine Is used whenever the program output is 
printed at times that do not coincide with the built-in program clock. 
Linear interpolation between two values (XI, Yl) and (X2, Y2) is used 
as illustrated in Figure 2-7. For a value X between XI and X2, the 
corresponding value Y is given by the linear relation: 



<Y 2 - v 


(52) 


The variable X represents the time and Y the variable to be interpolated. 


2.12 


VIEW (CONFIGURATION) FACTORS 


To take advantage of the axisymmetric nature of a solar 
receiver, a special view (or configuration, geometric) factor model is 
used. The possible node shapes would be a plane disc, plane ring, 
cylindrical ring, conical ring, etc., as shown in Figure 2-8. 

Each axisymmetric lode surface is entered by the user in the 
cylindrical coordinate system bv four coordinates r^ r^ in the radial 
direction and Z^, Z 2 in the axial direction, respectively. The radial 
coordinate (r) will be measured from the receiver centerline and the 
axial coordinate Z can be measured from the cavity aperture for cavity 
receivers or from the bottom of the external receiver connecting with 
the engine section as shown in Figure 2-9. The selection of the Z 
datum is arbitrary. 

Each line ring as shown in Figure 2-10 will be defined by 
its radius (r) and the corresponding axial distance (Z); thus the sets 
(r^, Z^), (r^, Z 2 ) correspond to the two line ring coordinates that con- 
stitute each node. The set (r^, Zj) will always refer to the bottom 
ring and (r^, Z 2 ) will always refer to the top ring. Two cases will 
be discussed as follows: 
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2.12.1 


Calculation of Ir.terior-to- Interior Node View Factor 


For two axisymmetric nodes i, j, as shown in Figure 2-11, 
each node will be presented by four coordinates (r^, r^, Z^, Z 2 ). The 
view factors from bottom surfaces to top surfaces are computed as 
follows: 


F(i, j) - PU.Jj) - F (i,j 2 ) 

A(j.) A(jJ 

- lor F (J i* 1} - Mir F ( v *> 


A(j.) A(jJ 

y~ [F (j 1 , i 2 ) - F (j r i x )] - [FC^.i^-F^,^) ) 


A(i) 


or 


F (i,j) = 


( A(i 2 ) F (ij,^) - AUj) F (iyjj)) - (A(i 2 )F(i 2 ,j 2 ) - 


AdjjFdj.jj)) 


(53) 


where 


A(i 2 ) - 7T r (i 2 ). 


A(i) - Trtrd^ + r(i 2 )] ^[r(i 2 > - rd^J 2 + [Z(i 2 > - Zd^T 


and 


A(ij) » it r (i^) 


The above four view factors F(i 2 ,jj)> F(i 2 ,j 2 >, F(ij.,j 2 )» and FCij^jj), 
between the plane discs forming node boundaries should be calculated 
first before the view factor F(i,j) between nodes (i) and (j) can be 
calculated. The view factor between two plane discs, as shown in 
Figure 2-12, can be expressed by the closed analytical form: 
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Figure 2-11. Calculation of View Factor, Interior-to- Interior Nodes 






( 54 ) 



where r^ and are the radii of the higher and lower discs, respectively 
and X is given by 

2 2 

X - 1 + (r H /r L r + (s/r L ) (55) 

To facilitate the view factor calculations, it is necessary to define the 
radii vector (R) with a dimension of (2 * NN) , where (NN) is the number 
of axisymmetric nodes that exchange energy by radiation. Each node will 
be represented by two elements of the vector (R) . For example, the i*"* 1 
node will have a lower radius located at the (21-1)^ element of (R) , 
and a higher radius located at the (21)^ element of (R). 


[R] 



node 1 
node i 
node NN 


[Z] - 


►node 1 

i 

►node i 

J|node NN 


Similarly, it is necessary to define the axial vector (Z) with a 
dimension of (2*NN) to present the axial distances of the lower and 
higher rings that bound each node. The calculation of the four elementary 
disc-to-disc view factors, DDF(4), previously mentioned in Equation (S3) 
is done using the 4-element vectors X(4) , RH(4), RL(4) , S(4). Use is 
made of the view factor reciprocity relations and the view factor 
summation rule to complete the calculation process of the matrix F(i,j). 
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Equation (53) is applicable to an interior node surface 
looking into another interior node surface above it in the Z direction or 
an exterior node surface looking into another exterior node surface above 
it in the Z direction. The sign of F (i,j) in Equation (53) is reversed 
if node (J) is below node (i) . For an interior node surface looking 
into another exterior node surface or vice versa, the Flag matrix V (i,j) 
is added to distinguish between these two cases and separate the view 
factor calculations. 

2.12.2 Calculation of Interior-to-Exterior Node View Factor 

For this case, a rough set of calculations is followed. A 
modification of the view factor preserved in Reference 8 for two con- 
centric cylinders of the same finite length, illustrated in Figure 2-13, 
is made to approximate the case of two dislocated conic rings of unequal 
lengths. 


Although the view factor calculations are approximated in 
this case, they could be changed by the user and replaced by other 
expressions if they are more accurate. The procedure followed is out- 
lined below. 


For the two concentric cylinder nodes (i) and (j) sketched 
in Figure 2-13 where the node (i) is the inner surface node and the node 
(j) is the outer surface node, the view factor F (i,j) is given by: 
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Figure 2-13. Two Concentric Cylinders of Same Finite Length 


where 


R - r (i)/r (j) 

L - £/r(j) 

2 —2 

A » L* + R - 1 
2 —2 

B - \T - R* + 1 


(57) 


i 



To allow for conic rings that are not of equal length and 
also dislocated axially with respect to each other, as shown In Figure 
2-14, the view factor In Equation (56) was modified to use average radii 
for the nodes (1) and (J), and an average length for both nodes. Using 
vectors R (2*NN) and Z(2*NN), the equivalent radii and length can be 
approximated as: 
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NODE(j) 



Figure 2-14. Approximation of Two Conical Nodes in View 
Factor Calculation 


r (j) ■ INR * J .1) » equivalent radius of node (j) 

r (i) ■ OTR ■ ^ - ^ * equivalent radiuc r*f node (i) (58) 

i « E Q VL , p( 21) - Z(2i-l )j ,| ~Z(2j) - Z(2j-l)j _ equiva.^.c 

For dislocated concentric nodes, the cosine of the angle of dis • ocation 
(KOST) is used to modify the view factor F(i,j) given by Equation (56) 
and written as: 

(OTR - INR) 

KOST - . = *■ (59) 

V (XX) + (OTR- INR) 
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where XX is the dislocation or shift given by 


XX 


Z(2j) + Z(21-l) 
2 



( 60 ) 


Furthermore, to account for the small slope of each conical 
ring, the view factor between two dislocated cylinders is multiplied by 
the cosine of the total angle between the two conical surfaces: 
cos (</>£ + <p J where: 


Cos (<J>^ + 4>j) 


Cos 



R(2i) - R(2i-1) 
Z(2i) - Z(2i-1) 


-1 

+ tan 


R(2j-1) - R(2j) ) 

Z(2j) - Z(2j-1 ) j 


( 61 ) 


It should be noted that the above view factors calculated 
for interior- interior nodes or interior-exterior nodes may still need some 
minor alterations by the user to correct for partial shading by other 
nodes. The latter is sometimes inevitable in view factor calculations 
and the user should make several trials to make sure that the rules of 
summation and reciprocity of view factors are applied correctly and the 
errors in satisfying these rules are minimized. 
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SECTION 3 


EXPLANATION OF INPUT DATA AND OUTPUT FORMS 
3.1 INPUT DATA 

To facilitate the use of the program, a set of tables is 
provided in Appendix B, with a brief description, to explain how each 
input variable shall be entered. Table 3-1 lists alphabetically all 
the input data, whether they are in free format or not. More data 
explanation can be found in Appendix C. 

The program language is written in MBASIC. The structure 
of the program is divided into two files and is given in Appendix D. 
These two files are the [FIRSTPART] file, which includes all dimension- 
ing and initialization statements and a large space for entering input 
data, and the [LASTPART] file which contains the set of computational 
statements. Merging of the two files in the program execution is done 
using the command [MERGE] in MBASIC. The first file, [FIRSTPART], 
should be loaded by the user using the command [LOAD] to enter all data. 
Tele-computer terminals are chosen for their fast communication with 
the user. Statement numbers 450-2000 are used to enter data, as will 
be explained in the example in the next section. 

Vectors are entered either in an element-by-element form or in 
scattered forms. For example, if the vector [ESOL] has 5 elements such 
that 


ESOL (1) »0.2 
ESOL (2) * 0.0 
ESOL (3) “0.5 
ESOL (4) -1.0 
ESOL (5) “0.0 

either one of the following statements can be typed: 
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Table 3-1. Input Variables 


VARIABLE 

DESCRIPTION 


ADJ 

Conduction/convection identifier 


AMBNN 

Number of nodes representing ambient air 


ASOL 

Nodal surface area for solar/IR radiation exchange 

ASRF 

Nodal surface area for conduction/convection heat transfer 

CONDUC 

Flag for conduction heat transfer in problem 


CONVEC 

Flag for convection heat transfer in problem 


DEN 

Node density 


DTOUT 

Time interval for printing out results 


EIR 

Nodal infrared emittance 


ESOL 

Nodal emittance in the solar band 


F 

View factor matrix 


FLOW 

Flag for fluid flow in problem 


FNDIN 

Number of nodes that represent the fluid inlet 
receiver 

section of 

FNDOT 

Number of nodes that represent the fluid outlet 
of receiver 

section 

HCONV 

Convective heat transfer coefficient 


HRIN 

Time sequence vector for time-varying variables 


HROUT 

Starting time for printing out program results 


I FLUX 

Solar flux incident on each node 


KCOND 

Nodal thermal conductivity 


LNGTH 

Conduction path length between nodes 


MASFLO 

Flow direction identifier 


MINDTU 

Minimum time increment for zero-capacity nodes 


MFM 

Mass flow rate or multiplier 


MNIT 

Maximum number of iterations in Newton-Raphson 

method 

NDIN 

Number of data points in time-varying input 
of HRIN vector 

dimension 

NN 

Total number of nodes 


POWR 

Internal power generated in each node 


R 

Radial coordinates of nodes 


SAIRR 

Flag for solar/inf rared radiation 


SIGMA 

Stefan-Boltzman radiation constant 


SORSIN 

Source/sink node identifier 


SPHT 

Nodal specific heat 


SS 

Flag for steady state/ transient solutions 


SSTEM 

Temperature data of source/sink nodes 


TEM 

Nodal temperature 


V 

View factor flag 


VOL 

Nodal volume 


VPRNT 

Flag for printing view factor matrix 


Z 

Nodal axial coordinates 
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> 500 


ESOL (1) - 0.2, ESOL (3) - 0.5, ESOL (4) - 1.0 (scattered form) 

or 

> 500 REAL ESOL (5) / 0.2, 0.0, 0.5, 1.0/ 

or 

> 500 DIM ESOL (5) / 0.2, 0.0, 0.5, 1.0/ 

Remember that statements in MBASIC can take more than one Instruction 

but separated with commas as shown above. Default values are always 

equal to zero. 

Overwriting data is also possible in MBASIC, and could be 
used to correct for typing errors. For example: 

> 600 ESOL (1) - 0.9, ESOL (3) - 0.5 

> 700 ESOL (3) - 0.4 

is equivalent to 

> ESOL (l) t = 0.9, ESOL (3) = 0.4 

For vector or matrix elements that have identical values, the 
following shorthand notation could be used in MBASIC: 

> 800 A(l) , A(2) , A(3) - 1 

> 900 B(2, 3) , B(3,4) , B(3,10) - 0.5 


is equivalent to 


A(l) 

- 1 

B(2,3) 

- 0.5 

A(2) 

- 1 

B(3,4) 

■ 0.5 

A(3) 

- 1 

B(3,10) 

■ 0.5 


Also, the [FOR] modifier can be used to assign values for vector or 
matrix elements having identical values. For example: 
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> 1000 C ( I , 3) 


2 FOR I - 7 TO 12 


Is equivalent to 

C(7,3) - 2 , C(10,3) - 2 

C(8,3) - 2 , C(ll,3) - 2 

C(9,3) - 2 , C(12,3) - 2 

To input a large NXN matrix. It is preferred to input only 
the elements that are non-zero in the scattered form. For some 
symmetrical matrices, only half the matrix is needed as input, and the 
other half is completed by the program structure. The solved example 
given in the next section illustrates how the data are collected, 
entered, and used to calculate the performance of a cavity receiver. 
More information about using the MBAS1C programming language can be 
found in Reference 9. 

3.2 OUTPUT FORMS 

The program results in the following output: 

(1) Applicability of the summation rule of view factors 
The output of E F (i,j) for the i C ^ node should be 

j 

either 1.0 or zero. Zero is only applicable for 
nodes that are assumed not exchanging heat by 
radiation. The user should make the necessary 
changes in the view factors after the program is 
first run to make sure that the summation rule is 
applied to all radiation nodes whether they are 
partially shaded or not. 

(2) View-factor matrix 

This output is presented if the user selects the 
value 1 for [VPRNT], The view factors are printed in 
column form for each node. Non- zero view factors are 
only printed. 
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Nodal energy residue 

This Is applicable to steady state solutions only. 
Equilibrium nodes should show a very small energy 
residue (around 10” 4 - lO"' 9 Btu for example) at 
steady state. Source/slnk nodes should give a non- 
zero energy residue with a proper sign. Positive 
energy residue means that there Is a net energy 
accumulation to the node and vice versa. As a check 
on computational accuracy, using the first law of 
thermodynamics, the algebraic sum of energy residues 
to source/sink nodes must equal the sum of all 
external energy entering the receiver through solar 
flux or other Internal heat generation. 

Accumulated energy extracted 

This variable is summed over all the time steps 
starting from the Initial time of computations up 
to the printout time (t) . The accumulated energy 
extracted is computed for both quasi-steady state 
and transient solutions. The daily accumulated 
energy extracted when the solar receiver is modeled 
over one day is useful in determining the integrated 
effects of transients. 

Accumulated solar energy into receiver 

This variable is also summed over all the time 
steps, starting from the initial time of com- 
putations up to the printout time (t). The 
accumulated solar energy received is computed 
for both quasi-steady state and transient 
solutions. 
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(6) Receiver efficiency 

Tli is is the ratio of the accumulated energy extracted 
from the receiver up to time (t) to the accumulated 
solar energy into the receiver up to time (t). 

For steady state solutions, the time interval is 
taken as one hour. 

(7) Modal temperature profile 

The temperature distribution of the receiver nodes 
at time (t) is given in tabular form. This output 
is desirable in observing the location of high 
temperature spots, outlet fluid temperature, and the 
determination of isothermal contours. 

(8) Nodal minimum transient time, 

For each node, the minimum time increment allowed 
for stable transient computations is printed. Also, 
the time increment (MINDTU) specified by the user to 
distinguish between zero capacity nodes and equili- 
brium nodes is reprinted. 

(9) Indirect Output results 

In addition to the above direct program printouts, 
the program can be used, either through several runs 
or by minor code modifications, to give a complete 
performance analysis of a particular receiver in 
the following areas: 

(a) Sensitivity Analysis 

This will include the sensitivity of solar 
receiver performance to variations occuring 
in the major receiver parameters such as nodal 
emissivity in the IR band (EIR), solar band 
emissivity (ESOL) , fluid mass flow rate (MFM), 
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nodal solar flux (1FLUX), ambient temperature 
[SSTEM (1)] or the fluid inlet temperature, 

TEM (FNDIN) . 

The effects on system efficiency and output 
can be studied in each case. 

(b) Itemization of Heat Losses 

The total heat loss from the receiver can be 
divided into 3 main parts: 1) a radiation 

heat loss through the aperture (for cavity- 
shaped receivers) and from the external boundary 
surface; 2) a convection heat to the ambient air 
outside the receiver and/or inside the receiver; 
3) a conduction heat loss through the insulation 
layers or because of lateral exchange between 
neighboring nodes. 

TVo runs could be made to quantize the magnitude 
of each part lost. The first run assumes zero 
thermal conductivity of the insulation layers or 
other layers that the user sees as causing a 
significant heat loss to the surroundings. The 
difference between the heat losses from this 
first run and the total heat losses would give 
the conduction loss. The second run assumes 
zero convective heat transfer coefficient between 
all the nodes exposed to ambient air and the 
ambient air node. The difference between the 
heat losses from the second run and the total 
heat losses would give the convection loss. 

The radiation loss will then be the third part 
left after subtracting the conduction and con- 
vection parts. 
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(c) Parameterization Study 

This study may include the effects on receiver 
performance because of changes made in receiver 
geometry, Interior piping arrangements, different 
working fluids or operating conditions. Studies 
of full-day performance under varying solar flux, 
inlet fluid conditions and ambient temperatures 
could be made with results specific to the 
receiver under study. Also, system optimization 
taking into consideration the receiver and 
engine performance is another indirect output 
of the receiver program. 
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SECTION 4 


CASE STUDY OF A CAVITY-TYPE RECEIVER 


In this section, the input and output procedure used in the 
computer program is applied to an experimental receiver. A brief 
description of the receiver, input data and the output results are 
presented . 


4.1 RECEIVER DESCRIPTION 


The receiver understudy was built at JPL in 1977 as a part of 
an energy program aimed at studying solar-thermal electric power plants. 
The receiver is placed at the focal spot of a 2.9-m (9.5-ft) diameter 
paraboloid dish; both are located at the JPL Table Mountain facility. 

The receiver is referred to as the Table Mountain receiver for identifi- 
cation. 

The receiver body Is a cavity-type composed of four primary 
parts as shown in Figure 4-1, These parts are: 1) a heat exchanger 
section consisting of 34 pairs of ii-tubes made of Inconel 718. The 
tubes are brazed to a head plate which is also made of Inconel 718; 

2) a ceramic receiver countour that forms the inside surface; 3) an 
aperture piece with changeable aperture sizes; and 4) an insulation layer 
in the form of a stainless steel canister filled with calcium silicate. 

The receiver performance was tested using helium as the 
working fluid at various inlet temperatures and pressures. Only one set 
of operating conditions is presented here in the receiver simulation. 

The overall receiver is cylindrical, 34.8-cm (13.7-in) in 
diameter and 12.7-cm (5-in) high. The working receiver aperture is 
5.08-cm (2-in) in diameter. The theoretical focal size (sun spot) is 
1.59-cm (0.625-in) and the actual focal size is 2.03-cm (0.8-in). 
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Figure 4-1. Cavity-Type Receiver Understudy 





4.2 


PREPARATION OF INPUT DATA 


Tables 4-1 through 4-13 Include all the receiver data required 
to run the program. The tables are constructed following the instructions 
given for each parameter in Appendix C. The receiver was discretieed into 
32 nodes as shown in Figure 4-2. Nodes 1, 2, 3, and 4 represent the 
Interior ceramic surface of the cavity. Nodes 5 and 6 represent the outer 
insulation layer. Node 7 represents the head plate-supporting structure. 

To make use of the view factor subroutine for axlsymmetric nodes, the 
34 pairs of U-tubes were treated as two coaxial annulus channels. Nodes 
9, 10, 11, 12, 13, and 14 represent the annulus channel through which 
the fluid enters, while nodes 16, 17, 18, 19, 20, and 21 represent the 
annulus channel through which the fluid exits. Nodes 8 and 15 represent 
the tip of the U-shape connecting the two annulus channels. Nodes 22 
through 28 represent the fluid in its passage through the tubes. Node 29 
represents the state of the working fluid leaving the receiver and node 
31 represents the state at entrance. Since the receiver inlet temperature 
was kept constant during the receiver operation, regardless of the exit 
temperature variations or ambient conditions, node 31 was considered as a 
source/sink node. Node 30 represents the head plate through which the 
fluid paths are grooved, and node 32 represents the ambient air temperature 
around the outer receiver body as well as the interior cavity. The 
following subsections describe the input data. 

4.2.1 Data for Table 4-1 

The problem at hand includes all four modes of heat transfer, 
i.e., radiation, thermal convection, conduction, and fluid convection. 
Therefore, the variables SAIRR, CONVEC, CONDUC, and FLOW were assigned 
the value 1. For a steady state solution, SST is assigned the value 1, 
keeping in mind that the program is capable of also handling transient 
solutions. With a total of 32 nodes, two of which are source/slnk nodes, 
the rest of the variables in Table 4-1 are assigned. Since only the 
steady state case was considered at one particular hour of operation, 
one element of the (NDIN) vector was given. The maximum number of 
iterations (MNIT) was assumed as 6, based on past experience in iteration 
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Table 4-1 Scalar Variables 


No. 

Variable 

Value 

Type 

1 

SAIRR 

1 

Integer 

2 

CONVEC 

1 

Integer 

3 

CONDUC 

1 

Integer 

4 

FLOW 

1 

Integer 

5 

SST 

1 

Integer 

6 

VPRNT 

1 

Integer 

7 

NN 

32 

Integer 

8 

ss 

2 

Integer 

9 

KINDTU 

1/60 

Real 

10 

AMBNN 

32 

Integer 

11 

FNDIN 

31 

Integer 

12 

FNDOT 

29 

Integer 

13 

FTIME 

1.0 

Real 

14 

NDIN 

1 

Integer 

15 

HNIT 

6 

Integer 

16 

HROUT 

0.0 

Real 

17 

DTOUT 

1.0 

Real 


Table 4-2. SORSIN, TEM, EIR, ESOL & POWR 






convergence. The minimum time increment (MINDTU) , below which nodes are 
considered zero-capacity nodes, was entered as 1/60-hour and is only 
applicable when the user requests the transient solution, 

4.2.2 Data for Table 4-2 

Nodes 31 (fluid inlet section) and 32 (ambient air) were the 
only source/sink nodes as indicated in the SORSIN column. The initial 
temperature for all nodes, excluding node 31, was taken as 25°C (77°F). 

The inlet fluid temperature was kept fixed at 344.44°C (652°F) . The 
emissivity values in the infrared range (EIR) or the solar band (ESOL) 
w le best estimates. No attempt was made to measure any of these 
eirissivities. For the aperture, or ambient air node (node 32), perfect 
black body values must be entered, i.e., ESOL = EIR = 1. Since there is 
no internal heat generation in any node, the (POWR) vector is set to zero. 

4.2.3 Data for Table 4-3 

Nodal specific heat, density, and thermal conductivity were 
abstracted from the literature for the materials used in the receiver 

2 

construction. Helium properties were taken at a pressure of 206,73 N/cm 
(300 psia) and an average temperature of 373°C (703. 4°F) . The volume of 
each node is calculated from the receiver dimensions in Figure 4-2. The 
volume of nodes 31 and 32 (infinite capacity nodes) were set arbitrarily 
at 1. Although the volume of each node is not needed for a steady-state 
solution, it was given for reference or for use in a transient solution. 
The surface area exposed to either direct solar radiation or to infrared 
energy exchange (ASOL) was also computed from the receiver geometry. 
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Table 4-3. SPHT, DEN, KCOND, VOL & ASOL 


Node 

SPHT 

(Btu/lb°F) 

DEN « KCOND 

(lb/ft J ) (Btu/hr ft°F) 

VOL 

(ft 3 ) 


ASOL 

(ft 2 ) 

1 

0 . 

12 

511 

6. 

5 

0.00364 

0.017429 

2 

0 . 

316 

140 

0 . 

83 

0.00592 

0.080505 

3 

0 . 

316 

140 

0 . 

83 

0.00509 

0 . 

12344 

4 

0 . 

316 

140 

0 . 

83 

0.01233 

0 . 

21764 

5 

0 . 

2 

11 

0 . 

,04 

0.14199 


0 

6 

0 . 

2 

11 

0 . 

,04 

0.22998 


0 

7 


t 


f 


i 

0.02119 


0 

8 




1 



0.00008 

0.047305 

9 







0.00012 

0.0705 

10 

0 . 

12 

51 

1 

6. 

,5 

0.00017 

0 . 

,10104 

11 






! 

0.00015 

0 . 

08722 

12 







0.00013 

0 , 

,07503 

13 






1 

0.00015 

0 . 

,088802 

14 







0.00008 

0, 

,04909 

15 







0.00005 

0 . 

,02727 

16 







0.00006 

0, 

,03466 

17 







0.00011 

0.06877 

18 







0.00011 

0.064117 

19 







0.00009 

0.054522 

20 







0.00010 

0,058179 

21 







0.00006 

0. 

,03633 

22 

\ 

f 





0.00097 



23 

1 . 

'24 





0.00116 



24 



0.0961 

0 

.1457 

0.00069 



25 







0.00059 

0 , 

.0 

26 





! 


0.00041 



27 







0.00078 



28 







0.00070 



29 





\ 


0.00070 



30 

0 

.12 

511 

6 

.5 

0.01615 

0 

.021814 

31 

1.24 

0.0961 

0 

.1457 

1.0 


0 

32 

0 

.24 

0.075 

0.015 

1.0 

0.021814 
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A. 2. A 


Radiation View-Factors 


Tables A-A and A-5 were constructed simultaneously to complete 
the radiation view factors. Nodes 5, 6, 7, 22, 23, 2A, 25, 26, 27, 28 
and 29 were excluded entirely fr<m the radiation exchange and were left 
blank. The elements of the view factor matrix that have a flag equal to 
[3] in Table A-A, i.e. , V(i,j) = 3, were marked (X) in Table A-5 and 
were given as: 

F(2,3) » 0.1 

F(3,A) = 0.1 

F(3,9) = 0.AA8 

F(13 ,16) = 0.083 
F(13,17) = 0.563 
F(13,18) = 0.051 

The above view factors were estimated after a preliminary run of the 
program to allow for shading by other nodes and to satisfy the summation 
and reciprocity rules. This procedure may or may not be needed in most 
cases, but was introduced in this document to familiarize the user with 
this iteration process. 

A. 2. 5 Conduction and Convection Between Nodes 

Table A-6 shows the distinction between nodes exchanging heat 
by conduction or by convection. The surface area in contact (ASRF) was 
computed for these two modes of heat transfer and is marked by (X) in 
Table A-7. The elements of the ASRF matrix are given below for the 
example receiver. Additional conditions imposed for filling out the 
tables are given in Appendix C, The English system of units is used 
throughout this example, but can be readily replaced by the SI (metric) 
system if necessary. 
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Table 4-4. View Factor Flag [V] 
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ibSu 





■I 


j 
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Table 4-5. View Factor Matrix [F] i 


NODE 
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ASRF 

(1,2) - 

0.11454 

ASRF 

(5,6) = 

1.01450 

tl 

(1,5) - 

0.04363 

II 

(5,7) = 

0.33540 

tf 

(1,32) - 

0.11454 

II 

(5,32) = 

0.39540 

II 

(2,3) - 

0.05760 

If 

(6,32) - 

2.59836 

II 

(2,5) - 

0.09817 

11 

(7,30) - 

0.10123 

It 

(2,32) - 

0.07592 

It 

(7,32) - 

0.48380 

It 

(3,4) - 

0.04014 

It 

(8,25) = 

0.04729 

If 

(3,5) = 

0.15272 

If 

(8,32) - 

0.04729 

II 

(3,32) - 

0.12217 

It 

(9,24) = 

0.07150 

II 

(4,5) = 

0.20180 

II 

(9,32) = 

0.07150 

II 

(4,7) * 

0.08727 

II 

(10,23) - 

0.10079 

It 

(4,30) = 

0.08050 

II 

(10,32) = 

0.10079 

It 

(4,32) = 

0.02220 , 

11 

(12,22) - 

0.07200 

It 

(11,22) - 

0.08727 , 

11 

(13,23) = 

0.08857 

II 

(11,32) - 

0.08727 

• 1 

(14,24) = 

0.05236 

II 

(15,25) * 

0.02730 

If 

(16,26) - 

0.03436 

II 

(17,27)- 

0.06872 

It 

(18,28) - 

0.06152 

II 

(19,28) - 

0.05236 

It 

(20,27) - 

0.05803 

II 

(19,32) = 

0.05236 

It 

(20,32) = 

0.05803 

II 

(21,26) = 

0.03436 

It 

(30,32) - 

0.12960 

It 

(21,32) - 

0.03436 


— 



To complete the computations of the thermal conductance needed for 
the conduction heat transfer, the matrix (LNGTH) in Table 4-8 has to 
be completed in full. However, the LNGTH matrix is completed for those 
nodes showing ADJ*1. The elements of the non-symmetric (LNGTH) were 
computed (in ft) and were given as tabulated below. 
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Table 4-8. [LNGTH] Matrix 





■□I 
mni 


ini 


DnimBSSSnSSSSn! 


£»■■■■■■■■■■■! 


IS! 


IS! 








LGNTH 

(1,2) 

a 

0.0167 

LNGTH 

(2,1) 

m 

0.0375 

»t 

(1,5) 

- 

0.0583 

It 

(5,1) 

- 

0.0917 

ti 

(2,3) 

- 

0.0375 

It 

(3,2) 

at 

0.0583 

i« 

(2,5) 

St 

0.0417 

If 

(5,2) 

9E 

0.0917 

it 

(3,4) 

s 

0.0583 

ft 

(4,3) 

as 

0.1125 

ti 

(3,5) 

a 

0.0208 

It 

(5,3) 

- 

0.0917 

it 

(4,5) 

as 

0.0260 

ft 

(5,4) 

= 

0.0917 

ti 

(4,7) 

= 

0.0333 

If 

(7,4) 

* 

0.0917 

n 

(5,6) 


0.0917 

1! 

(6,5) 

at 

0.0917 

it 

(5,7) 

* 

0.2083 

ft 

(7,5) 

= 

0.0333 

ti 

(4,30) 


0.1104 

ft 

(30,4) 

= 

0.0375 

it 

(7,30) 

= 

0.0333 

ft 

(30,7) 

= 

0.03755 


For thermal convection, the convective heat transfer coefficient (HCONV) 

was assembled as shown in Table 4-9. The heat transfer coefficient was 

2 

taken as approximately 1 Btu/hr°F ft between node 32, which is the 

ambient air, and nodes 1, 2, 3, 4, 8, 9, 10, 11, 19, 20, and 21. For 

nodes 5, 6, 7, and 30, the convective coefficient was considered as 

2 

approximately 2 Btu/hrF ft to allow for the unaccounted radiation loss 
from the receiver exterior surface. Note that the above values could be 
better refined if more analytical information was given about the con- 
vection mechanism inside and outside the receiver cavity. For the con- 
vective coefficient between the working fluid (helium) and the receiver 
tubes, exemplified in the set of nodes (8,25), (9,24), (10,23), (11,22), 
(12,22), (13,23), (14,24), (15,25), (16,26), (17,27), (18,28), (19,28), 
(20,27) and (21,26), the following formula for turbulent flow inside 
tubes was used: 

\7 no- 1 „ 0.8 _ 0.4 

Nu * 0.023 R Pr 

n 
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Table 4-9. [HCONV] Matrix 


Hgggggggggg^^ 


!■■■■■■■■■■■■■■■■■■■■ 

• -BBBBBBBBBBBB* BBBB ■ 

^^^■bBBBBBBBBBBBBBBBBBBIBBBBB^H 
HMKflBBBBBBBBBBBBBBBBBnBBBBBB^H 
rai ifiiBBBBBBBBBBBBBBBBBDllllill^M 

^■bbbbbbbbbbabbbbbbbbbdbbbbbbb^H 
HbbbbbbbbbbbbbbbbbbbEbbbbbbbbIIH 

^■BBBBBBBBBBBBBBflBBBBBBBEBBBBBBBBl 


IBBBBBI 

IBBKJBBI 


ILBI 

in 

in 


in 


where Nu is the Nusselt number, Pr is the Prandtl number and Rn is the 
Reynolds number. The equivalent, or "Hydraulic" diameter (D^) to be used 
in the above formula is determined from the expression 

D • 4x cro8S sectional area 
h wetted perimeter 

For an annulus area of an inside diameter d. and an outside diameter d , 

i o 

the hydraulic diameter is found to be (d -d.). 

o 1 

In the present example, helium gas flows at a rate of 42.11 kg/hr 

2 

(92.9 lb /hr) at 206.73 N/cm (300 psia) pressure. Since the inlet helium 
temperature was 344. 4°C (652°F) , the property of helium was taken arbi- 
trarily at an average temperature of approximately 372. 8°C (703°F) based 
on about 55.6°C (100°F) temperature rise through the receiver tubes. The 
Reynolds and the Prandtl numbers were 52,400 and 0.659, respectively. 

The convective heat transfer coefficient was then calculated as 
2 

580 Btu/hr°F ft . For different helium flow rates, the convective 
coefficient (H) , in Btu/hr°Fft^, is written in the form 

H - 15.453 m f °* 8 

where m^ is the flow rate in lb/hr inside the receiver tubes and H is the 
convective coefficient between the helium gas and the tubes in the 
example receiver. 

4.2.6 Fluid Convection 

The flow direction matrix (MASFLO) was developed as in 
Table 4-10. The closed fluid loop (only one in this example) was 
considered to go from nodes 31-22-23-24-25-26-27-28-29-31. Either zero 
or 1 is entered to identify the flow direction. In Table 4-11, 
variations of the mass flow rate (MFM) with time (HRIN) were allowed. 

Only one value of the flow rate was imposed in the fluid loop (92.9 lb/hr 
in Table 4-11) and was assumed to remain constant thereafter. Table 4-11 
was completed by entering the fixed temperature for source/sink nodes 
SSTEM(l) and SSTEM(2) . Note the two sources/sink nodes in this example 
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TV: 


I 


Table 4-10. [MAS FLO] Matrix 


Node 



were the ambient temperature (node 32) and the inlet fluid (node 31) , 
respectively. 


Table 4-11. SSTEM(l) , SSTEM(2) and MFM 


HRIN 

■ 

2 

3 

■ 

5 

6 

l 

■ 

8 

9 

10 

11 

SSTEM 

(1) 

77 

1 

1 

1 

1 

1 

1 

1 

1 



SSTEM 

(2) 

652 

1 

| 

1 

1 

1 

1 

1 

1 



MFM 

92.9 

1 

1 

1 

1 

1 

1 

1 

1 




4.2.7 Flux Distribution 

The solar flux distribution i.. 'able 4-12 was predetermined 
using an optics computer program developed for the on-going solar -thermal 
power system analyses. The intensity concentration ratio (c) divided by 
the parabolic mirror reflectivity (p) is given in Figure 4-3 together 
with a sketch of the cavity receiver made to fit the optics program. 

Table 4-13 was used to determine the solar flux on nodes ?, 8, 19, 20, 

21, and 30 in terms of the values given in Figure 4-3. 

The total solar energy retained inside the cavity receiver was 

found from a previous calorimeter test to range from 4.097 kW^ 

(13,98L>.7 Btu/hr) to 4.47 kW t (15,256.1 Btu/hr) . Arbitrarily taking the 

lower end of the scale, the equivalent solar intensity falling on the 

2 

mirror surface would be 13,983.7/83.07 or 168.3 Btu/hr ft . Note that 
the latter is always less than the actual solar intensity falling on 
the mirror to allow for the receiver intercept factor caused by the 
energy spill-over from the limited receiver aperture size. The last 
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Table 4-12. Solar Flux [IFLUX] 
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Figure 4-3. Solar Flux Distribution on Example Receiver 
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Table 4-13. Solar Flux Distribution 


Node 

C/p 

c (1) 

asol (2) 

ft 2 

C*AS0L 

iflux (3) 

Btu/hr ft 2 

3 

300 

270 

0.12344 

33.33 

45,441 

8 

800 

720 

0.047305 

34.06 

121,176 

19 

0 

0 

0.054522 

0.0 

0.0 

20 

100 

90 

0.058179 

5.24 

15,147 

21 

200 

180 

0.03633 

6.51 

30,294 

30 

200 

180 

0.021814 

3.93 

30,294 

Total 




83.07 



(1) Assuming a mirror reflectivity (p) of 0.9 

(2) Taken from Table 4-3 

2 

(3) Product of (C) column multiplied by 168.3 Btu/hr ft 


column In Table 4-13 Is thus completed and the data entered back into 
Table 4-12. Again, in Table 4.12, only one column is used corresponding 
to one set of solar flux values which remain constant thereafter. The 
user has the option of completing Table 4-12 for a transient solution by 
repeating this procedure at preselected time intervals during the solar 
day. 


4.2.8 Nodal Coordinates 

In Table 4-14, the axial and radial coordinates of each node 
were given based on a reference point taken arbitrarily at the cavity 
aperture center with the Z-direction taken along the receiver canterline. 
Each node in the cavity interior is represented by two bounding circular 
rings. Therefore, four coordinates (r^, r^, Z^, and Z are needed. 

Only axisymmetric nodes that exchange heat by radiation need to be con- 
sidered to make use of the special view-factor subroutine embedded in 
the code. For example, nodes 5, 6, 7, 22-29, and 31 were assumed not 
exchanging heat by radiation to neighboring nodes and their coordinates 
were entered as zero. 

After completing the input data tables, (LOAD) the [FIRSTPART] 
to enter the data, then (MERGE) the [LASTPART] of the program and the 
program is ready for execution using the (RUN) command. 

4.3 ANALYSIS OF PROGRAM RESULTS 

A copy of the program results is shown in Figures 4-4 and 4-5. 
The output for the steady state solution starts by a list of the radiation 
view factors (if requested by the user), a check on the summation rule of 
the radiation view factors, the accumulated thermal energy extracted by 
the receiver, the accumulated solar energy incident on the receiver 
interior, the receiver efficiency, and the nodal residual energy. The 
user may repeat this exercise using different values for the receiver 
parameters for the performance evaluation, the transient solution and 
the sensitivity analysis. 
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Table 4-14, Nodal Coordinates 


Radial Axial 

Coordinates Coordinates 


Node 


QH 


r 2 


Z 1 


Z 2 

1 

0, 

.0833 

0.0833 

0 


0, 

.0333 

2 

0 , 

.0833 

0, 

,1583 

0 , 

.0333 

0 , 

.1083 

3 

0 . 

.1583 

0 . 

,175 

0 , 

.1083 

0 , 

.225 

4 

5 

0 , 

.175 

0, 

,1333 

0 , 

1 

.225 

| 

0 , 

1 

.4458 

6 

0 

0 

0 , 

.0 

0 , 

.0 

0 , 

.0 

7 



j 

\ 

J 

i 



8 

0, 

.14167 

0, 

,0708 

l 

0 , 

.1458 

0 

1458 

9 

0, 

.14167 

0, 

.14167 

0 , 

.1458 

0 , 

.225 

10 

0 , 

.14167 

0.1333 

0 , 

.225 

0 , 

.341667 

11 

0. 

,1333 

0 . 

.1333 

0 , 

.341667 

0, 

.4458 

12 

0. 

.11667 

0 . 

,1125 

0.341667 

0, 

.4458 

13 

0. 

.125 

0. 

,11667 

0 , 

.225 

0 , 

.341667 

14 

0 . 

,125 

0, 

,125 

0 , 

.1625 

0 , 

,225 

15 

0. 

,0833 

0, 

.125 

0 , 

.1625 

0, 

.225 

16 

0 . 

,0833 

0 . 

,091667 

0, 

.1625 

0, 

.225 

17 

0.091667 

0. 

,095833 

0 , 

.225 

0.341667 

18 

0. 

,095833 

0 . 

.10 

0, 

.341667 

0 , 

,4458 

19 

0. 

,08333 

0. 

.08333 

0 , 

.341667 

0, 

.4458 

20 

0 . 

.075 

0. 

,08333 

0, 

.225 

0, 

.341667 

21 

22 

23 

24 

0 , 

.0708 

0 . 

,075 

0 , 

.1458 

0 , 

.225 

25 

26 

27 

28 
29 

0 , 

,0 

0 . 

.0 

0 , 

.0 

0, 

.0 

30 

0.0833 

0 . 

,0 

0, 

.4458 

0 , 

.4458 

31 

0 , 

.0 

0 , 

.0 

0 , 

.0 

0 , 

.0 

32 

0.0 

0, 

.0833 

0, 

.0 

0 , 

.0 
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ENTER ZERO FDR TRANSIENT ANALYSIS DR ENTER 1 FDR QUASI STEADY STATE ONE*! 

ENTER NUNBER DF NODES 1 32 

ENTER NUNBER DF SOURCE-SINK NODES *2 

ENTER ND OF TINE DEPENDENT DATA POINTS* I 

ENTER ZERO IF THERE IS NO RADIATION EXCHANGE ELSE ENTER 1*1 

ENTER ZERO IF THERE IS NO CONDUCTION ELSE ENTER 1*1 

ENTER ZERO IF THERE IS NO CONVECTION ELSE ENTER 1*1 

ENTER ZERO IF THERE IS NO FLUID FLOW ELSE ENTER 1*1 

ENTER 1 TO PRINT VIEW FACTORS ELSE ENTER 0*1 
NON-ZERO ELENENTS OF THE VIEW FACTOR MATRIX ARE * 


F < 

1 B 

1> = 

. 1303526 

F < 

1 9 

2) = 

.23487491 

FC 

1 9 

3> = 

• 8501 0817E-01 

F < 

1 9 

8) * 

. 69947446E-01 

F< 

1 9 

19) = 

• 69776335E-02 

F < 

1 9 

20) = 

• 1924146E-01 

F < 

1 9 

21) = 

• 43644969E-01 

F < 

1 9 

32) = 

.40995018 

F < 

2b 

1> = 

.50349029E-01 

F (. 

2b 

2> * 

.36236416 

F < 

^9 

3> = 

.1 

F < 

2b 

9> = 

. 33987289 

F < 

2b 

19) = 

.26434368E-01 

F < 

2b 

20) = 

. 401 18394E-01 

F< 

2b 

21) = 

. 20979597E-01 

F < 

2b 

32) = 

.93815902E-02 

F < 

3 b 

1) = 

. 12002878E-01 

F < 

3b 

2) = 

• 6521 7652E-01 

F < 

3b 

3> = 

.12858836 

F < 

3 b 

4) = 

. 1 

F < 

3b 

9) = 

.448 

F < 

3b 

10) = 

.15747265 

F < 

3 b 

11) = 

. 38718466E-01 

F < 

4b 

3) = 

• 56718564E-01 

F< 

4b 

4) = 

. 36997634 

F < 

4b 

9)* 

.6955166E-01 

F < 

4b 

10> = 

.24778143 

F < 

4b 

11) = 

.25597202 


F < 

8b 

1)* 

• 25770986E-01 

F < 

8b 

2)* 

.66349273 

F< 

3b 

8) = 

. 17794836 

F< 

3b 

32) = 

.13278794 

F < 

9b 

3) = 

.78442779 

F < 

9b 

4) * 

.21471245 

j- < 

9b 

9) * 

• 8597672E-03 

F < 

10b 

3) * 

. 19238226 

F < 

10b 

4) « 

.53370762 

F < 

10b 

1 0> ■ 

.27391014 


Figure 4-4. Program Output for Steady-State Solution 
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P< 11, 3 >* 

PC 11, 4) = 
PC 11, U>* 

PC 12, 12) . 
PC 12, 15) = 
PC 12, 16)* 
PC 12, 17)« 
PC 12, 13) a 


. 12556642 
.63374333 
.23563371 

.16354773E-01 

• 361 32731E-01 
• 7661 7332E-01 
.12373056 

• 7416146 


PC 13, 13>» 
PC 13, 15) a 
PC 13, 16> = 
PC 13, 17) a 
PC 13, 13) a 

PC 14, 14) a 
PC 14, 15)a 
PC 14, 16)= 
PC 14, 1 7) a 
PC 14, 13)= 

P < 15, 12) = 
PC 15, 1 3> = 
PC 15, 14)= 
PC 15, 15>= 
PC 15, 16)= 
PC 15, 17) = 

PC 16, 12)= 
PC 16, 13) = 
PC 16, 14>= 

P C 16, 15>= 
PC 16, 16) = 


.21002237 
. 32377627E-01 
. 33000001E-01 
.563 
• 51E-01 

. 46035633E-01 
.27035323 

• 42027565 

• 1337507 
•73513676E-01 

. 33407533E-01 

• 30274573 
•43661 377 

. 7*061 472E-01 
. 32320032E-01 
.2351 3606E-02 

. 16535463 
•21264722 
.53513723 
• 25323343E-01 
• 47701537E-03 


PC 17, 12) = 
PC 17, 13) = 
F C 17, 14)= 
PC 17, 15>= 
p C 17, 17)= 

PC 13, 12)= 
PC 13, 13) = 
PC 13, 14)= 
PC 13, 13)= 

PC 13, 1> = 
PC 13, 2) = 
PC 13, 13>= 
PC 13, 20)= 
PC 13, 2D* 
PC 13, 30> = 
PC 13, 32)= 

PC 20, 1)= 

FC 20, 2)« 

P C 20, 13) = 
PC 20, 20) = 

P< 20, 21 >* 

P C 20, 30) * 
PC 20, 32) = 


■1405136 
.72703333 
■ 13116652 
• 3 325422 3E- 03 
. 34335533E-03 

.36736137 
•70635563E-01 
. 6037376E-01 
•62270463E-03 

.22305277E-02 
• 330321 33E-01 
•41705336 
.22273707 
. 3374546 3E-01 
.27635037 
. 32433323E-02 

•57642363E-02 
. 5551 3343E-01 
•20373574 
■43126326 
.15453312 
•30413273E-01 
. 1 371 1542E-01 


PC 21, 

1)» 

PC 21, 

2>* 

PC 21, 

13>« 

PC 21, 

20) « 

PC 21, 

21 >« 

PC 21, 

30> « 

PC 21, 

32) * 


. 2033324E-01 
. 46431357E-01 
•50645703E-01 
.24757053 
•56364613 
. 26023033E-01 
. 33677247E-01 


Figure 4-4. 


Program Output for Steady-State Solution 


(Cont'd) 


4-27 


F < 

30* 

19>* 

F< 

30* 

20> * 

F< 

30* 

21>* 

P< 

30* 

30>* 

F< 

30* 

32)* 


F < 

32* 

1>* 

F < 

32* 

2> * 

F< 

32* 

3> * 

F< 

32* 

19>* 

F < 

32* 

20> * 

FC 

32* 

21> * 

F< 

32* 

30>* 

F < 

32* 

32)* 


.21462495 
•43377136E-01 
•1664716E-01 
. 32670736E-01 


. 32776329 
. 34646522E-01 
.23315663 
. 20613733E-01 
• 36594155E-01 
• 66121 704E-01 
• 32670735E-01 
.1934232 


CHECK ON 

SUMMATION RULE OF F 

1 


. 99999999 

0 


1 

1 


1 

1 


1 

0 


0 

0 


1 

TIME <HR> 

HCC 

ENSY EXTRCBTU) 

0 


0 

NODE NO 

NODE 

TEMP* F 

1 

77 


2 

77 


3 

77 


4 

77 


5 

77 


6 

77 


7 

77 


3 

77 


9 

77 


10 

77 


11 

77 


12 

77 


13 

77 


14 

77 


15 

77 


16 

77 


17 

77 


13 

77 


19 

77 


20 

77 


21 

77 


22 

77 


23 

77 


24 

77 


25 

77 


26 

77 


27 

77 


23 

77 


29 

77 


30 

77 


31 

652 


32 

77 



1 

1 

1 

1 

0 


* cc 


SDL EN6Y 
0 


1 

1 

1 

0 

0 

INCBTU) 


0 

1 

1 

0 

0 

RECVR EFF 


0 

1 

1 

0 

1 


Figure 4-4. Program Output for Steady-State Solution (Cont’d) 


4-20 


NODE NQ 

NODE RESIDUE* BTUH 

1 

-.19579337E-01 

2 

-.3331 1336E-01 

3 

-.44436155 

4 

.14346313 

5 

33146973E-05 

6 

0 

? 

• 12207031E-03 

3 

. 65429638E-01 

3 

.16366577 

10 

. 72337695E-01 

11 

. 43413334E-01 

15 

. 23693307E-03 

13 

50232473E-03 

14 

.47350603E-03 

15 

-. 1 3664365E-03 

16 

- • 2354 0 062E - 0 3 

1? 

-.21 362305E-03 

13 

-. 37976456E-04 

13 

. 7 33233 3 3E- 02 

20 

. 1 1 192322E-01 

21 

• 70037342E-02 

22 

. 1 353125E-02 

23 

0 

24 

1 3531 25E-02 

25 

0 

26 

• 13531 25E- 02 

2? 

0 

23 

0 

23 

0 

30 

• 1625061E-02 

31 

11006.521 

32 

2377. 1373 


T I ME (HR > HCC ENSY EXTR(BT'J> HOC SDL ENSY IN(BT'J> RECVR EFP 

1 1 1 006. 521 13?33.»73 .73703303 


NODE NQ 

NODE TEMP* F 

1 

1331.1556 

5 

1413.2723 

3 

1426.3722 

4 

346.57669 

5 

301.31039 

6 

36. 160637 

p* 

747.33626 

3 

331.5369 

3 

742.96092 

10 

634. 05167 

11 

673.32035 

12 

660.37732 

13 

672.17061 

14 

690.90236 

15 

727.10236 

16 

735.37341 

17 

743.53195 

13 

746. 99315 

13 

761.43443 

20 

773.15372 

21 

734.41602 

22 

660.3262 

23 

671.7915 

24 

690.67439 

25 

727.35152 

26 

735.71692 

27 

744.05597 

23 

747.54605 

23 

747.54605 

30 

735.71452 

31 

652 

32 

77 


> 


Figure 4-4. Program Output for Steady-State Solution (Cont'd) 
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>RUN 

ENTER ZERO FOR TRANSIENT ANALYSIS OR ENTER 1 FOR QUASI STEADY STATE ONE» 0 

ENTER NUMBER V NODES* 32 

ENTER NUMBER OF SOURCE-SINK NODES* S 

ENTER NO OF TINE DEPENDENT DATA POINTS* 1 

ENTER ZERO IF THERE IS NO RADIATION EXCHANGE ELSE ENTER 1*1 

ENTER ZERO IF THERE IS NO CONDUCTION ELSE ENTER 1*1 

ENTER ZERO IF THERE IS NO CONVECTION ELSE ENTER 1*1 

ENTER ZERO IF THERE IS NO FLUID FLOW ELSE ENTER 1*1 

ENTER 1 TO PRINT VIEW FACTORS ELSE ENTER 0*0 
CHECK OH SUMMATION RULE OF F 


1 

.99999999 1 

1 

0 

0 

1 1 

1 

l 

1 

1 1 

1 

l 

1 

1 1 

0 

0 

0 

0 0 

0 

0 

0 

1 



NODE nd 

NODAL WIN DT»MR 



1 

.8766495E-01 



2 

• 85866 027E-01 



3 

.23618223 



4 

.19601015 



5 

.24215351 



6 

. 933856 03E- 01 



7 

. 1 0889373 



8 

. 17551506E-03 



9 

. 17326226E-03 



10 

.17621335E-03 



11 

.17553484E-03 



12 

. 18341 398E-03 



13 

. 17652056E-03 



14 

.1651 081 36-03 



15 

. 17690475E-03 



16 

.177915996-03 



17 

.176148256-03 



18 

. 1 8354969E-03 



18 

• 18303406E-03 



20 

.176304596-03 



21 

.186544116-03 



22 

.554151516-06 



23 

.616595286-06 



24 

.437653736-06 



25 

.44433766-06 



26 

.311328016-06 



27 

.494062436-06 



28 

.462664816-06 



29 

.727943226-06 



30 

.976057216-01 



31 

.598637526-06 



32 

. 12675796-05 




Figure 4-5. Program Output 

for Transient 

Solution 


4-30 


HIM TINE INCREMENT IN SECONDS IS 60 

TIME<HR> ACC ENGY EXTR <BTU> ACC SQL ENGY IN<BTl» RECVR EFF 
0 0 0 0 

NODE NO NODE TENP»F 


1 

77 

2 

77 

3 

77 

4 

77 

5 

77 

6 

77 

7 

77 

3 

77 

9 

77 

10 

77 

11 

77 

12 

77 

13 

77 

14 

77 

13 

77 

16 

77 

17 

77 

18 

77 

19 

77 

20 

77 

21 

77 

22 

77 

23 

77 

24 

77 

23 

77 

26 

77 

27 

77 

28 

77 

29 

77 

30 

77 

31 

632 

32 

77 


Figure 4-5. Program Output for Transient Solution (Cont’d) 
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TIME <HR> ACC ENSY EXTRCBTU) 
. 3E-01 463.5103 


ACC SDL ENSY IN<BTU> 
633. 13364 


RECVR EFF 
.66333141 


NODE NO 

NODE TEMPfF 

1 

116.45512 

2 

219.65307 

3 

334.96397 

4 

131.55163 

5 

73.642746 

6 

77. 00396 

7 

30.461258 

3 

356.15227 

9 

719.16234 

10 

631.11226 

11 

669.99344 

12 

658.06944 

13 

665.94315 

14 

630.04375 

15 

713.50502 

16 

721.45219 

1? 

729.14673 

13 

732.10132 

19 

744.54955 

20 

756.75453 

21 

763.49033 

22 

657.60997 

23 

665.61737 

24 

679.33221 

25 

713.71633 

26 

721.75096 

27 

729.55602 

28 

732.56934 

29 

732.56934 

30 

97.47731 

31 

652 

32 

77 


Figure 4-5. Program Output for Transient Solution (Cont'd) 
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TIME <HR> ACC ENGY EXTR <BTU> 
.1 989. 5876? 


ACC SDL ENGY IN<BTU> 
1399.3673 


NODE ND 

NODE TEMP#F 

1 

197.88085 

8 

343.63176 

3 

646.58393 

4 

186.80766 

5 

34.679968 

6 

77.063084 

7 

91.803678 

3 

356.38851 

9 

780.43358 

10 

681.68148 

11 

670.89976 

18 

653.16585 

13 

666. 17894 

14 

630.55818 

15 

714.05444 

16 

788. 00546 

17 

789.70544 

18 

738.66374 

19 

745.13544 

80 

757.34004 

81 

769. 06984 

88 

657.70885 

83 

665.35153 

84 

630.34053 

85 

714.86765 

86 

788.3058 

87 

730.11714 

88 

733.13487 

89 

733.13487 

30 

113.58588 

31 

658 

38 

77 


RECVR EFF 
.66476647 


Figure 4-5. Program Output for Transient Solution (Cont'd) 
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APPENDIX A 


PROGRAM FLOW CHART 


A-l 














1 




3 











o 


COMPUTE ACCUMULATED 
ENERGY AND EFFICIENCY 

ENEX 

SEIN 

EFF 



1 NO 

( *• ) 


A- 6 


GOSUB 

INTERPOLATION 

SUBROUTINE 



INTERPOLATION SUBROUTINE 






































APPENDIX B 
LIST OF VARIABLES 


B-l 


NAME 


DESCRIPTION 


DIMENSION 


UNITS 


A 

View factor parameter 

— 

— 

*ADJ 

ADJ (I,J) is conduction/convection 

identifier between i^ and j** 1 
nodes such that 

ADJ(I,J) * 1 means that nodes i,j, 
are exchanging heat 
by conduction, 

ADJ(I,J) * 2 means that nodes i,j 
are exchanging heat by 
convection 

(NNxNN) 

matrix 


*AMBNN 

The number of nodes representing 
ambient air 

— 

— 

*ASOL 

th 

ASOL (i) is the i node surface 
area between nodes i,j for 
radiation heat transfer only. 

(NN) 

vector 

ft 2 

*ASRF 

ASRF(i,j) is the surface area between 
nodes i,j for conduction/convection 
heat transfer only, 

■ 0 diagonally 

■ 0 for 2 non-ad jacent nodes 

ASRF is only given if the problem 
contains conduction or convection. 

(NNxNN) 

matrix 

ft 2 

B 

View factor parameter 

— 

— 

CAPV 

CAPV(i) is the thermal capacitance 
for the i t ^ 1 node. 

(NN) 

vector 

Btu/°F 

CLOCK 

Main program time counter. Clock 
is incremented after each interval. 

Scalar 

hr 

CLOCK1 

Time counter for interpolation sub- 
routine . 

Scalar 

hr 

CONDC 

CONDC(i,j) is the conductance 
matrix between nodes i,j for 
conduction heat transfer 

(NNxNN) 

Btu/hr°F 

*CONDUC 

Flag for conduction heat transfer 
in solution 

— 

— 


*Varlables that must be given by the user as input data. 


B-2 


NAME 


DESCRIPTION 


DIMENSION 


UNITS 


CONVC 

CONVC (i,j) is the conductance for 
nodes i,j exchanging heat by 
convection 

(NNxNN) 

Btu/hr°F 

*C0NVEC 

Flag for convection heat transfer in 
solution 

— 

— 


CONVEC * 0 No convection 

■ 1 Solution includes con- 
vection 



DDF 

Disc-to-disc view vector for axi- 
syrranetric configuration factor 
subroutine 

(4) 

vector 

— 

*DEN 

DEN(i) is the density of the i^ 
node material 

(NN) 

vector 

lb/ft 3 

DTAU 

Minimum time step (At) for transient 
solution. Also used as the time step 
for quasi-steady state solution. 

Scalar 

hr 

DTAUV 

DTAUV (i) is the minimum time step of 

the i node, calculated from 

stability subroutine. Applicable 
to transient solution only. 

(NN) 

vector 

hr 

DTEM 

DTEM(i) is the i ^ nout new 
temperature at (t+At) 

(NN) 

vector 

°R 

*DTOUT 

Time interval to print out results. 

Scalar 

hr 

EFF 

Accumulated receiver efficiency 
from starting time till time (t). 

— 

— 

*EIR 

EIR(i) is the infra-red emittanca 
of the i t ^ 1 node surface. 

(NN) 

vector 


ENEX 

Accumulated energy extracted from 
the receiver, from starting time 
till time (a) 

Scalar 

Btu 

ENEXIN 

Accumulated energy extracted from 
the receiver, used for interpolation 
subroutine. 

Scalar 

Btu 

ENEX1 

Accumulated energy extracted from 
the receiver at time (t^) 

Scalar 

Btu 


Variables that must be given by the user as input data. 


B-3 


NAME 


DESCRIPTION 


DIMENSION 


UNITS 


EQVL 

Equivalent length for view 
factor calculation. 

Scalar 

ft 

*ESOL 

ESOL(i) is the emittance of the 
th 

i node in the solar radiation 

band 

(NN) 

vector 


El, E2, E3 
E4, E5, E6 

Parameters for view factor sub- 
routine at v = 2. 

— 

— 

*F 

F(i,j) is the view (or config- 
uration) factor between nodes 
i,j that "see” each other. 

(NNxNN) 

matrix 

— 

FCF 

FCF(i,j) is total view factor 
matrix (F) between nodes i,j 
for the infrared radiation 
exchange . 

(NNxNN) 

matrix 


FCFST 

FCFST(i,j) total view factor 
matrix (£*) for solar radiation 
exchange between nodes i,j. 

(NNxNN) 

matrix 

" 

FLO 

FLO(i,j) is the fluid conductance 
between nodes i,j. 

(NNxNN) 

matrix 

Btu/hr°F 

*FLOW 

Flag for fluid flow in solution 

Flow =0 No fluid flow 

= 1 Fluid energy is 
considered in 
sol ution 



FLUX 

FLUX(i) , direct solar flux 
incident on i tn node. 

(NN) 

vector 

Btu/hr ft 

*FNDIN 

Number of node which represents 
the fluid inlet section of 
receiver 

— 


*FNDOT 

Number of node which represents 
the fluid outlet section of 
receiver. 


“ “ 

*FTIME 

Final clock or run time 

Scalar 

hr 


^Variables that must be given by the user as input data. 


t-T- 


B-4 


NAME 


DESCRIPTION 


DIMENSION 


UNITS 


FU 

FU(i) is the sum of view factor 
from node (i) to all other nodes 
in enclosure = ZF(i,j) 

j 

Transfer matrix between excitation 
and response vectors for IR 
radiation 

(NN) 

vector 


FI 

(NNxNN) 

matrix 

— 

FIS 

Transfer matrix between excitation 
and response vectors for solar 
radiation exchange 

(NNxNN) 

matrix 

— 

HEAP 

Heat Energy Analysis Program title 

— 

— 

*HCONV 

HCONV(i,j) is the convective heat 
transfer coefficient between 
adjacent nodes i,j exchanging 
heat by convection 

(NNxNN) 

matrix 

Btu/hr ft 

*HRIN 

HRIN(k) is the time in hr corres- 
ponding to the k t ^ 1 input transient 
data , 

(NDIN) 

hr 

*HROUT 

Starting time for printing out 
program results 

Scalar 

hr 

1,11,11 

Index 

— 

— 

IDN 

Identity matrix (Kronecker 
delta of size NN) 

(NNxNN) 

matrix 

— 

I FLUX 

IFLUX(i,k) is the input solar 

flux incident on the 

t h th 

i node at the k time interval 

(NNxNDIN) 

Btu/hr ft 

INR 

Inner radius for view factor 
calculation 

Scalar 

ft 

INTER 

Flag for interpolation subroutine 

Scalar 

ft 

IREXC 

IREXC (i,j) is the infra-red 
radiation exchange coefficient 
between nodes i,j 

(NNxNN) 

matrix 

CM 

LW 

IREXC1 

Special matrix formed from IREXC 

(NNxNN) 

matrix 

ft 2 


^Variables that must be given by the user as input data. 
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NAME 


DESCRIPTION 


DIMENSION 


UNITS 


JAC 

K, KK 
*KCOND 

KD 

KOST 

L, LL 
LIREXC 

*LNGTH 

LR 

*MASFLO 


*Variable 


Index 


Jacobian matrix for steady state 
solution iteration 


JAC(iJ) 


3 RES(i) 
3 TEM(j) 


Index 

KCOND(i) , conductivity of 
node material 

Kroencker delta (identity matrix) 

Parameter for view factor sub- 
routine 


Index 

Linearized matrix coefficient 
IREXC to suit the stability 
subroutine 

LNGTH(i,j) is the distance 
between nodes i,j exchanging 
heat by conduction measured 
from center of node i to 
separating surface between i.j 

Parameter for view factor sub- 
routine 

MASFLO(i,j) is the fLow direction 
identifier between nodes L,j 

- 0 No fluid exchange 

= 1 Flow is from node j to 

node i for a single 
fluid loop 

* Mass flow rate from node ' 
to node i for more than 
one fluid loop 


(NNxNN) 

matrix 


(NN) 

vector 

(NNxNN) 

matrix 


(NNxNN) 

matrix 


(NNxNN) 

matrix 


Scalar 


(NNxNN) 
r.at r ix 


that must be given bv the usv i as input: data. 


Btu/°F 


Btu/hr ft°F 


Btu/hr°F 


ft 


ft 


NAME 

DESCRIPTION 

DIMENSION 

UNITS 

MINDTU 

Minimum time interval below which 
nodes are considered zero-capacity. 
For transient solutions. 

Scalar 

hr 

*MFM 

The Mass flow rate (assuming one 

fluid circuit in the receiver) at 

, , th _ 

the k time interval, lor more 

than one fluid loop, MFM is taken 

as 1. 

(NDIN) 

vector 

lb /hr 

*MNIT 

Maximum number of iterations in 
Newton-Raphson method. 

— 

— 

N 

Index 

— 

— 

*NDIN 

Number of data points in time 
varying input or dimension of 
HRIN vector. 

— 

— 

NE 

Number of equilibrium nodes 
NE = NN-SS 

— 

— 

*NN 

Total number of nodes 

— 

— 

OTR 

Parameter for view factor sub- 
routine 

Scalar 

ft 

PH 1 1 

Parameter for view factor sub- 
routine. 

Scalar 

radian 

PHI J 

Parameter for view factor sub- 
routine . 

Scalar 

radian 

PI 

The number tt = 3.1416 

— 

— 

*P0WR 

POVR(i), internal heat generated in 

i th node, by electrical resistance, 
nuclear reaction, etc. 

(NN) 

vector 

Btu/hr 

*R 

Radial coordinates for node 
surfaces exchanging radiation. 

(2*NN) 

vector 

ft 

RES 

RES(i), residual of i ^ node, or 

t h 

sum of heat input to i node 

(for quasi-steady state solution) 

(NN) 

vector 

Btu/hr 


*Variables that must be given by the user as input data. 
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NAME 


DESCRIPTION 


DIMENSION 


UNITS 


RH 

Radius of higher disc (or ring) 
forming the surface of the axi- 
symmetric node in DDF view factor. 

(4) 

vector 

ft 

RHO 

RHO(i) , spectral reflectivity of 

i^ node in the infrared radiation 
spectrum 

(NN) 

vector 


RHOST 

RHOST(i), spectral reflectivity of 

i node in the solar radiation 
band 

(NN) 

vector 

mm 

RL 

Radius of lower disc (or ring) 
forming the surface of the axi- 
symmetric node in DDF view factor. 

(4) 

vector 

ft 

RR 

Parameter for view factor 
subroutine. 

Scalar 

ft 

*S 

Height between higher disc and 
lower disc in DDF view factor. 

(4) 

vector 

ft 

*SAIRR 

Flag for solar and infrared 
radiation existence in solution. 

— 

— 


SAIRR « 0 no radiation 
= 1 solar & IR 

exchange are present . 



SEIN 

Accumulated solar energy incident 
on receiver from start to time (t). 

Scalar 

Btu 

SEININ 

Accumulated solar energy incident 
on receiver used for interpolation 
subroutine. 

Scalar 

Btu 

SEIN1 

Accumulated solar energy incident 
on receiver from start up to (x^). 

Scalar 

Btu 

*SIGMA 

Stefan-Boltzman radiation constant 
1.714 X 10" 9 

Scalar 

Btu/hr ft 2e4 R) 

SOLCON 

SOLCON(i), net solar radiation 
exchange by i t ' 1 node. 

(NN) 

vector 

Btu/hr 


*SORSIN Source/sink node identifier (NN) 

SORSIN(i) « 0 equilibrium node vector 

* 1 source/sink node 


Variables that must be given by the user as input data. 
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NAME 


DESCRIPTION 


DIMENSION 


UNITS 


*SPHT 

*SS 

SSK 

**SST 

**SSTEM 

T 

*TEM 

TEMIN 

TEMNEW 

TEM1 

UNM 

UNV 

V 


^Variables 


SPHT(i), specific heat of i^ node 
material. 

(NN) 

vector 

Btu/lb 

Number of source /sink nodes 

— 

— 

Steady state index to stop 
temperature calculation when 
input data no longer changes. 


— 

Flag for transient and steady state 
solutions identifier. 

— 

— 

SST = 1 steady state 
= 0 transient 



Source/sink temperature profile 

(SSxNDIN) 

matrix 

°R 

Index 

— 

— 

TEM(i) is the temperature of i^ 
node 

(NN) 

vector 

°R 

Nodal temperature vector used for 
interpolation subroutine 

(NN) 

vector 

°R 

TEMNEW(i), updated temperature for 
i node at (t 4* At) . 

(NN) 

vector 

°R 

Nodal temperature vector at end 
of time (t^). 

(NN) 

vector 

°R 

Matrix with all elements = 1 

(NNxNN) 

matrix 

— 

Vector with all elements = 1 

(NN) 

vector 

— 

V(i,j) is a radiation exchange 
flag between nodes i,j 

(NNxNN) 

— 

V(i,j) == 0 nodes i,j do not "see" each other 

= 1 interior cone surface (i) to 

interior cone surface (j) 

= 2 interior cone surface (1) to 



exterior cone surface (j) 

= 3 asymmetric node F values to 


be entered by the user 
* 4 nodes i, j to be calculated 
using summation rule 

that must be given by the user as input data. 
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NAME 


DESCRIPTION 


DIMENSION 


UNITS 


*VOL 

VOL(i), volume 

of i C ^ node 

(NN) 

vector 

ft 

*VPRNT 

Flag for view factor matrix 
printing 

— 

— 

X 

Parameter vectc 
calculation 

*-*♦ [I] 

>r for DDF 

(4) 

vector 


XX 

Parameter for radiation subroutine 

Scalar 

ft 

*z 

Axial coordinates of lower and 
higher discs forming node (for 
axisymmetric view factor subroutine) 

(2XNN) 

vector 

ft 

ZCN 

Zero-capacity node flag. ZCN is 1 
for a zero-capacity node; else 
ZCN(i) is 0. 

(NN) 

vector 


ZCN1 

Zero-capacity nodes identifier in 
transient solution. 

— 

— 


*Variables that must be given by the user as input data. 
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APPENDIX C 
INPUT DATA TABLES 


C-l 


SAIRR 


Flag for solar-infrared radiation heat transfer 

SAIRR ■ 0 if problem does not include nodes 
exchanging solar & infrared 
radiation 

else SAIRR ■ 1 


CONVEC 


Flag for conduction heat transfer 

CONDUC «0 if problem does not include nodes 
exchanging heat by conduction 

else CONDUC - 1 


CONDUC 


Flag for thermal convection heat transfer 

CONVEC * 0 if problem does not include nodes 
exchanging heat by convection 

else CONVEC - 1 


FLOW 


Flag for flow convection 

FLOW - 0 if problem does not include 
fluid nodes flowing from one 
node to another 

else FLOW - 1 
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I SST 

I I 

Flag for transient/8teady state solutions 

SST “ 1 For steady state solution 
SST ■ 0 For transient solution 


AMBNN 


Number of node representing the ambient air 

VPRNT 



Flag for printing the view factor matrix 

VPRNT ■ 1 the program will print the 
view factor matrix 


else VPRNT ■ 0, No printing 


FTIME 


Final clock or run time of program 
simulation, in hours 


MINDTU 


Minimum time interval ( At ) in hours 
below which nodes are considered sero- 
capaclty nodes. MINDTU is only required 
in transient solutions (SST ■ 0). 10/3600 
could be used as a first trial. 










NN 


Total number of nodes In the problem including 
source/sink nodes and equilibrium nodes 


NDIN 


Number of data elements entered in HRIN table, 
or dimension of HRIN vector 


ss 


Number of source/sink nodes in problem. These 
nodes can gain or lose any amount of heat without 
affecting their temperature 


MNIT 


Maximum number of Newton-Raphson Iterations 
required. Usually not less than 6 is recommended. 
Accuracy of computations may not be affected very 
much with increasing MNIT. 
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HROUT 


Starting time (in hours) at which program results 
are printed 


dtout 


Time interval (in hours) for printing program 
results 


FNDIN 


Number of node representing the fluid inlet 
section. If no fluid circuit is present:, enter 
zero 


FNDOT 


Number of node representing the fluid outlet 
section. If no fluid circuit is present, enter 
zero 
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Table C-l. SORSIN 


NODE 

_1_ 

_2_ 

El 

□ 

a 

a 

□ 

a 

□ 

ID 

D 

m 

m 

14 

IS 

m 

IQ 

ns 

B3 

20 

SORSIN 














1 

1 

1 

1 

l 

1 

1 


Source/ 8 ink node identifier. SORSIN shall be either zero or 1 
according to the following: 

SORSIN ■ 0 For an equilibrium mode 

SORSIN • 1 For a source/sink node 

Nodes that represent ambient air, constant temperature or a phase 
changing process are examples. 


Table C-2. TEH 


NODE 


_2_ 

3_ 

4_ 

_5_ 


7 

8 

_9_ 

10 

11 

12 

13 

m 

IS 

m 

IS 

ID 

□ 


i 














1 

1 

l 

l 

1 

1 

1 


Initial temperature nodes in °F. The temperature vector is 
used as an Initial estimate for Iteration in steady-state 
solutions or as an initial value for transient solutions. 
Default value is 77*F for all the nodes. If the initial 
temperature of certain node(s) are different from 77*F, the 
user shall enter the new value (s). 







Table C-3. EIR 


NODE 

K 

K 

□ 

a 

1 

□ 

B 

a 

□ 

□ 

m 

12 

ts 

14 


m 

m 

m 


20 1 

UJ 

1 

1 

1 








1 








1 



The semi-grey emissivity of each node 'n the infrared region 
extending approximately from a wavelength of 4 microns (4000 A) 
to infinity. For nodes where the emissivity is unknown, too 
small, or does not exchange heat by radiation, a value of zero 
shall be entered. 


Table C-4. ESOL 


NODE 

a 

□ 

a 

a 

B 

B 

B 

B 

□ 

IQ 

0 

m 

is 

D 

IS 

m 

m 


Q 

20 

ESOL 

l 

l 

l 

l 

1 

1 

1 

1 

1 

1 

1 

l 

l 

1 

1 

l 

l 

1 

1 



The semi-grey emissivity of each node in the short wave (solar) 
region below 4 microns (4000 A) wavelength. For nodes where 
the emissivity is unknown, too small, or does not exchange heat 
by radiation, a value of zero shall be entered. 










Table C-5 


POWR 


NODE 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

POWR 






















Internal power generated in each node in Btu/hr, due to 
chemical reaction, nuclear excitation or heat generated by 
electrical resistances. Default value is zero for each 
element. 


Table C-6 . SPHT 


NODE 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

SPHT 






















The specific heat of each node material is in Btu/lb°F. Default 
value is zero for each element. 
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Table C-7. DEN 


NODE 

D 


□ 

D 

□ 

□ 

□ 

□ 

B 

m 

D 

12 

m 

14 

IS 

m 

□ 

ES 

IS 

20 

DEN 

1 


1 

1 

1 

1 

1 

1 

1 

l 

1 

1 

1 


1 

1 

1 

1 

1 

1 


The density of each node material in lb/ft^. DEN is only 
needed for transient solutions. For steady state solutions, 
DEN table shall be left blank. Default value is zero for 
each element. 


Table C-8. KCOND 


NODE 

n 

□ 

El 

a 

H 

n 

B 

B 

B 

m 

ID 

12 

m 

m 

IS 

m 

IS 


m 

20 

KCOND 

l 

l 

1 

l 

1 

l 

1 

1 






l 

l 

l 

l 

l 

l 

1 


The thermal conductivity of node material, in Btu/hr ft°F. 
Default value is zero for all nodes. 
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Table C-9. VOL 


| NODE 

I 

2 

El 

0 

a 

□ 

□ 

□ 

□ 

D3 

D 


i 


IS 

m 

□ 

m 

m 

20 

VOL 



1 

1 

1 

1 

1 

1 

1 

1 

1 




1 

1 

1 

l 

l 



The volume of each node material in ft^. VOL is needed for 
transient solutions only. Default value for each element is 
zero. Source/sink nodes shall be assigned an arbitrary 
volume, e.g., 1 ft^. 


Table C-10. ASOL 


NODE 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

m 

20 

ASOL 



















1 



Surface area in ft^ of a node subject to heat transfer by 
radiation with other neighboring nodes, or a node receiving 
direct solar radiation. The ASOL will be internally computed 
by HEAP when the axisyimnetric view factor subroutine is used. 
Otherwise, the user shall enter the values of ASOL. For nodes 
that do not exchange heat by radiation completely, a zero 
value shall be entered. 











Table C-ll. View Factor Flag [V] 


I 


i 


NODE 



View factor flag {V] is an integer to identify the viewing configuration 
between any two nodes. Only half of the [V] matrix (hatched area) shall 
be entered to include the diagonal elements. The flag element V (i f j) 
represents how the i^ node "sees" the j** 1 node as follows: 


V (i, j) » 0 

V (i, j) - 1 

V (i.j) - 2 

V <i,j) - -2 

V (i, j ) - 3 

V (i,j) • 4 


Nodes (i) and (j) do not "see*' each other either directly 
or because of shading or obstruction by other nodes. 

Node (i) is the interior of a conical ring, exchanging 
radiation with node (j), which is the interior of a 
conical ring. 

Node (i) is the interior of a conical ring, exchanging 
radiation with node (j), which is the exterior of a 
conical ring. 

Node (k) is the exterior of a conical ring exchanging 
radiation with node (j) which is the interior of a 
conical ring. 

Nodes (i) and (j) do not fit an axisymmet ric model 
described above in V - 0, 1, 2, or -2 and have a given 
view factor. 

Users free choice node-node where view factor will be 
determined by the summation rule. Each row in the 
hatched [V] half matrix shall contain one node only 
with V « 4. 
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View factor F (i,j) between nodes (i) and (j) is only entered whenever 
the flag V (i,j) is 3; the case of two nodes where the view factor is 
not determined by HEAP subroutines. Only elements in the hatched area 
need to be entered. Non-applicable cases shall be left blank. 
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Table C-13. ADJ 


NODE 



ADJ is a "flag" matrix to identify the type of heat transfer between 
adjacent and non-ad jacent nodes. ADJ is an integer that can be as 
follows : 

Blank or 0 For non -ad jacent nodes 

1 For two adjacent nodes exchanging heat by conduction 

2 For two adjacent nodes exchanging heat by convection 

Only half of the matrix needs to be entered (the hatched area) because of 
symmetry. Diagonal elements shall be left blank. Default value is 
zero for each element. 


C-13 







































Table C-14. ASRF 



Contact surface area between two adjacent nodes (in ft^) that exchange 
heat by conduction (ADJ * 1) or convection (ADJ = 2), Default value 
is zero* Non-applicable boxes shall be left blank. Only half of the 
matrix (shaded area) shall be entered with blank diagonal elements. 
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Table C-15 . LNGTH 


NODE 



D 

2 

3 

Q 

5 

6 

D 

8 

9 

10 

a 

12 

13 

14 

15 

a 

17 

18 

19 

a 

21 

22 

23 

24 

25 

mm 


























Li_ 




















































wm 


























5 


























6 


























mm 


























8 


























9 


























10 


























11 


























12 












■ 














13 












i 














14 

















■ 

■ 








15 




, 






















16 


























17 





















■ 





18 





















i 





19 





















■ 





20 


























21 


























22 


























23 


























24 


























25 



























The heat transfer average path length between two adjacent nodes 
exchanging heat by conduction. LNGTH Is entered in ft for nodes with 
ADJ * 1 only. LNGTH (i,j) is measured from the center of node (i) to 
separating surface between nodes (i) and (j). 


C-15 











































T.iMc i!- 1 €» . 


IICONV 


i 

! 

I 



NODE 



Coefficient of heat transfer by convection in Btu/hr ft°F. HCONV is 
entered for nodes exchanging heat by convection only (ADJ >2). Only 
half of the matrix, excluding the diagonal elements, shall be entered. 
Non-applicable boxes shall be left blank. Default value is zero for 
all elements. 
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Table C-17. MAS FLO 


NODE 



3DS3S3SEESEESEBSBSESS 




: 


MASFLO (i,j) is an Identifier for the direction of flow between nodes 
(1) and (j) if a single fluid loop is present In the receiver. MASFLO 
(i,j) is also an identifier for both the direction and mass flow rate 
if the receiver includes more than one fluid loop. For a single fluid 
loop, MASFLO (i,j) is entered as 1 if the flow direction is from node 
(j) to node (1) and zero if the flow direction is from node (1) to node 
(j). For a single fluid loop, MFM (Table C-18) will be the mass flow 
rate in lb/hr. For receivers with more than one fluid loop, the actual 
mass flow rate in lb/hr will be entered in MASFLO (i,j) if the fluid 
flow direction is from node (j) to node (1), and zero otherwise. In 
the latter case, MFM shall be entered as 1. 


i 


C-17 



































Table C-18. HRIN, SSTEM and MFM 


HRIN 

D 

2 

3 

0 

5 

D 

B 

0 

B 

10 

11 

12 

13 

B 

15 

16 

S 

1 

19 

S 

21 

22 

S 

24 

SSTEM 

(1) 

1 

1 

l 

1 

l 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

SSTEM 

(2) 

1 

1 

1 

1 

1 

1 

1 

1 

fl 

fl 

fl 

B 

B 

1 

fl 

B 

fl 

B 

fl 

1 

1 

1 

1 

1 

MFM 

1 

1 

1 

1 

1 

1 

1 

1 

B 

B 

1 

B 

B 

1 

B 

1 

fl 

1 

1 

1 

1 

1 

1 

1 


HRIN Tine domain data, giving the sequence of clock tines (not 

necessarily of equal intervals) when some variables change 
or are not to be re-specified by the user. 

SSTEM Source/slnk temperature variations with time domain, entered 
in °F. SSTEM (1) gives the temperature of the ambient node, 
and SSTEM (2) gives the temperature of the inlet-to-receiver 
fluid node. Both the ambient node and inlet -t,»-receiver node 
are considered source/sink nodes here, but could be changed 
in the program structure. 

MFM Mass flow multiplier in lb/hr. This is the mass flow rate of 

the working fluid inside the single receiver loop. Variations 
of MFM with time are entered according to the number of data 
points in the time domain. For receivers with more than one 
loop, the number 1 should be entered as a flow baseline and 
variations of flow rates with time should be entered as the 
ratio of new flow rate to the base line. Tables C-17 and 
C-18 should be completed simultaneously. For receivers with 
more than one fluid loop, the mass flow rate of each loop is 
assumed to vary simultaneously with time so that their ratios 
remain unchanged. 
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Table C-19. IFLUX 


NODE 


9 

D 

2 

3 

4 

5 

Q 

D 

D 

D 

S 

5 

12 

13 

14 

15 

16 

17 

s 

19 

s 

21 

22 

23 

s 

25 

mm 













■ 

■ 











r 

2 













i 

■ 












3 


























mm 


























L!j 


























mm 


























mm 


























n 


























KB 


























10 


























11 

■ 

■ 

■ 

■ 

■ 

■ 

■ 

■ 

■ 

■ 

■ 

■ 

■ 

■ 

■ 

■ 


■ 

■ 

■ 

■ 

■ 

■ 

■ 


12 


















■ 

■ 

■ 

■ 

■ 

■ 

■ 


13 


























KOI 


























15 


























16 


























17 


























18 


























19 



























20 

















— J 









21 















—I 


“i 

i 









22 























_ 



23 


























24 















"" 




— J 







25 



























Solar radiation intensity in Btu/hr ft^. Variations of IFLUX with time 
intervals are provided in the tables. The values from 1 to 24 in the 
HRIN row are time intervals. 
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Table C-20. Nodal Coordinates 



RADIAL 

COORDINATES 

AXIAL 

COORDINATES 

NODE 

r 1 

r 2 

2 i 

*2 

1 





2 





3 





4 





5 





6 





7 





8 





9 





10 





11 





12 





13 





14 





15 





16 





17 





18 





19 


t 



20 






This nodal coordinates table Is needed only when the view factor matrix 
Is calculated using the axlsymmetric node subroutine. Each axlsymmetric 
node (whether It is a conical or cylindrical ring, or a planar disc) is 
defined by four coordinates: (r^,Z^) for the bottom boundary, and (r 2 ,Z 2 ) 
for the top boundary where Z 2 > Z^. The radial distance r is measured 
from an arbitrary datum. For nodes that are not exchanging heat by 
radiation, the coordinates shall be left blank. Default value is zero 
for each coordinate. All values are entered in ft. 


SURFACE NORMAL 



C-20 


APPENDIX D 
PROGRAM LISTING 


D-l 


FILE 1 

[ FIRSTPART ] 


D-2 


100 INPUT USING CHAPC13)+' ENTER ZERO FOP TRANSIENT ANALYSIS 
OR ENTER 1 FOR OUASI STEADY STATE DNE»»'«SST 
ISO INPUT USING CHAR C 1 3> ♦ ' ENTER NUMBER OF NODES I O ' I NN 
140 INPUT USING CHAR fl3> ♦ ' ENTER NUMBER OF SOURCE-SINK NODESm'lSS 
160 INPUT USING CHAR Cl 3> ♦'ENTER NO OF TIME DEPENDENT DATA POINTS««' JNDIN 
180 INPUT USING CHAR C13) ♦'ENTER ZERO IF THERE IS NO RADIATION EXCHANGE 
ELSE ENTER lits'sSAIPR 

S00 INPUT USING CHAR C 1 3> ♦ ' ENTER ZERO IF THERE IS NO CONDUCTION 
ELSE ENTER lSw'JCONDUC 

SSO INPUT USING CHAP C 13) + ' ENTER ZERO IF THFRE IS NO CONVECTION 
ELSE ENTER IStt'iCONVEC 

S40 INPUT USING CHAR Cl 3) ♦' ENTER ZERO IF THERE IS NO FLUID FLOW 
ELSE ENTER 1 !«' I FLOW s 

260 INPUT USING CHAP Cl 3) ♦' ENTER 1 TO PRINT VIEW FACTORS ELSE ENTER 0s»'sVPRNT 
SSO CLOCK=0« SI6.MA=1 . 714E-09 

300 f DIMENSION AND INITIALIZE ALL ARRAYS 

320 REAL ASPF CNN.NN) .FIR CNN) .ESOL CNN) » F CNN* NN) < FCF CNN. NN) . FCFST CNN. NN) . 

RHD CNN) . RHOST CNN) . IREXC CNN. NN> . IREXC1 CNN. NN) . L IREXCCNN. NN) , 

POWR CNN) » CONDC CNN. NN) . HCONV CNN. NN) . UNM CNN. NN) =1 . (CD CNN. NN) s I DN CNN) r, 
LNGTHCNN.NN) .KCONDcNN) 

340 REAL CONVCCNN.NN) .FI CNN.NN) » F 1 S CNN. NN) . TEM CNN) . TEMNEW CNN) . DTEM CNN) , 

TEM1 CNN) , TEM IN CNN) . MASFLO CNN.NN) .FI O CNN, NN) , CAPV CNN) , SPHT CNN) , 

VOL CNN) . DEN CNN) , AD J CNN , NN) , SOL CON CNN) , ASOL CNN) . DTAUV CNN) , 

I FLUX CNN « ND I N ) .SSTEMCSS.NDIN) .FLUX CNN) .HRINCNDIN) , MFM CNDIN) 

360 REAL X C4> » RH C4) , RL C4 > , S C4) » DDF C4) , R <2»NN) . ? <3*NN) . FU CNN) , JAC CNN. NN) , 

V CNN. NN) . UNV CNN) =1 , SOPS IN CNN) . RES CNN) , ZCN CNN) 

380 TEM Cl) =77 FOP 1=1 TO NN 

400 ! THE FOLLOWING SPACE IS FOR ENTERING DATA 

CHECK THF PROGRAM MANUAL 


420 HPOUT=O.DTOIJT=1 , FNP I N=3 1 . AMBNN=32 » FNDOT =29 

440 MN I T=6 . FT I ME= 1 

460 MINDTU* 1-^60 

430 S0RSINC31) » SDPSINC32) =1 

500 TEM <3 1) =652 

520 REAL EIRC32) '.75, .75. .75. .75,0. .5.. 5. .5. .5. .5, .5,. 5. .5. .5.. 5, .5,. 5, .5. 
. 5* . 0« 0? 0> 0» 0« 0» .5» 1/ 

540 PERL ESOL . 156* . 156* . 156* . 156* Or .5* .“it .5> .5» .5» ,5» .5* 

, 5 * ■ 5 i ,5**5* • 5 • , 5 r Hi On fli fl> On 0 ? Or Or , 5 * Mr 1 •' 

560 PERL tPHTO$>/. 1 £r .316* . 316b .316. .£* .£. . 1£. . 1£. • 1£, .1£* . 1£* . 1£* 

. 1 £ * . 1 £ <• . 1 £ * . 1 £ * . 1 £ r . 1 c r * 1 £ * , t £ * . 1 £ * 1 . £4 . ! . £ 4 * 1 • £4 * 1 . £4 * 1. £4 * 

1 . £4r 1 . £4* 1 . £4. . 1£* 1 . £4 r . £4 

590 RFRL DFN* ?£> -'51 1 r 1 4 0, 140r 140* 1 1 * 1 1 r 51 1 r 51 1 . 511 * 5 1 i • 51! * 51 1 * 51 1 * 

51 1 r 51 1 .51 1 *SJ t r 51 l r5U r5U *51 1 * . 0961 r . 0961 * . 0961 r . 0961 * . 0961 r 
. 0961 * . 096 1 . . 0961 *511*. 0961 • . 075' 

600 PERL KCOND * 3 £ > ’'6. 5* . 83* . $3* . $ 3 « . 04r • 04 1 6, 5. 6 • 5*6. 5r 6. 5* 6. 5* 6. 5 r 
6. 5*6. 5 *6. 5 *6. 5*6. 5*6. 5* 6. 5*6. 5 *6,5* . 1457r • 1457r • 1457* . 1 457 r 
. 1457b . 1457 r . 1457 r . 1457r6.5r . 1457r . 015 
6£ 0 V *1 R £ '• • V *■ l R 3 > R V *• 1 B 9 > R V *■ t R 1 9 > R V *■ 1 R £ O') * V •' 1 * £ 1 > = 1 
►,4 0 V • 1 r 3£ ) r V •’ t- * 9 ) r V <i £ r 1 9 •* • V ■■ c’ r £ 0 * * V •' £ r £ 1 1 = 1 
660 V<£* 3£> . V- Xr » rV'IPr 15> r V< 1 ,3* 15 > r V *1 4 r 1S> =;1 

6*5 0 V '• 1 9 • £ 0 ) r V *' 1 9 * £ 1 ) r V • 1 9 « 3 fi > « V t 1 9 * 93 > r V •' £ 0 • £ 1 1 r V *• £ 0 * 3 0 * ► V * £ Q r 3 £ 1 — 1 

700 V • £ 1 . ?.0 > r V • £ 1 r 3£ > t V <"i 0. 3$ > = 1 

7£0 V*: r 1 in , v •:?* 1 l » • V* 4 r9:. , V <4r 1 O ' r V* 4r 1 1 :• -£ 

74 0 V < 1 £ 9 16)* V • 1 £ r 1 7> * V •: 1£* 19) * V* >14* 16) * VI 4* 1 7) * V* 14r 1 *’ £ 

760 V '15* 16) * V* 15* 1 7> -£ 

79 0 V 1 • 1 « * V »: £ • £ > * V •' 3 » 3 ) • V •' 19*19)* V •• £ 0 r £ 0 ) * V * £ 1 * £ 1 ) * V < 3£ r 3£ > r *■,•' • 8 * 8 • • 

V •• 9r 9 ) r V 1 Or 1 0 • r V ♦ 1 1 * t 1 > r V < 1 £ * l £ > r V • 19* 1 9 > r V <14* 1 4 • * V • 15,1 5 « =4 
$00 V •: 1 6 r 16> r V < 17* 1 7> * V 19* 1 9) *4* V* 4* 4> * V • ^ m r 30 » =4 
8£0 V4* 1 0 > =£r V *. 30* 3£>=1 

94 m V 1 3 * 4 ) * V * £ * ? > » 3 • P 1 3 * 4 » * . 1 * P ♦ £ * 3 1 “ • 1 
$6 0 •: 3 r 9 > * 3 < F 1 »<n = . 44 .-: 

$80 V< 1 8 r Ik) r V •: 13*1 7) * V *13* 1 3 > »9tF* 1 ?r 16' = . 083* F •: 13* 17) *. 563* F • 1 3* 1$> *. 051 
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900 REAL R<64> 03 33* . 03 33* . 0333* . 1533* . 1583. . 175* . 175* . 1 333* 0* 0* 0* 0* 0* 0* 

. 14167. . 0703. . 14167. . 14167. . 14167. .1333* . 1333. . 1333. . 1 1667. . 1 123. 

. 123* . 1 1667. . 123* . 123* . 03333. . 123. . 0333. . 091667. . 091667. . 095333. 

. 093333* . 1 * . 03333* . 03333. . 073. . 08333* . 0703. . 075. 0. 0. 0. 0. 0. 0. 0. 0. 
0*0»0*0*0»Q. 0.0*. 03 33.0. 0*0.0.. 083 3-< 

920 REAL 2 <64> ✓ O*. 0333.. 0333.. 1093*. 1083.. 223.. 223. . 4453. 0*0. 0*0. 0*0. 

. 1453* . 1438* . 1459. .225* .225. .341667. .341667, .4438. .341667. 

.4453* .225* . 341667, . 1625, .225, . 1625, . 1 625 ». 1625, .225, .225, 
.341667* . 341667* .4458* .341667, .4458* .225, .341667* .1458, .225, 

0*0 *0*0*0, 0*0, 0*0, 0*0*0, 0*0, 0,0,. 4459 * . 4458^ 

940 ADJC1 ,2) *ADJ<1 ,5) *ADJ<2»3) , ADJ<2,5> *ADJ<3,4) , ADJ (3*5) ,ADJ<4»5)*1 
960 MU <4* 7> » AD J <4* 30) , ADJ (5*6) , ADJ (5*7) *ADJ <7, 30) *1 < 

930 ADJ U* 3?> =2 FOR 1=1 TO 21 

1000 AD J <9. 25> , ADJ <9, 24 > *ADJ<1 0,23) , ADJ Cl 1 »22> , ADJ <12, 22) ,ADJ<13*23>=2 
1020 ADJ<14*24> ,ADJ<15,25> *ADJ<16*26) , ADJ<17*27) ,ADJ<18,29) *ADJ<19,28>=2 
1040 ADJ <20* 27> , ADJ <2 1 * 26> *ADJ (29*30) ,ADJ(30*31> , ADJ (30* 32) =2 
1 06 0 AD J (29* 30) , AD J ( 3 0 * 3 1 ) * 0 
1080 ADJ(I* 32> =0 FOR 1 = 12 TO 19 

1100 ASRF ( I » 2) » ASRF < 1 * 32) “ . 1 1454* ASRF (1 *5) =. 04363* ASRF (2,3) =. 0576, ASRF < 

2,5>=. 09917, ASRF <2, 32>*. 07592, ASRF (3,4) =. 04014, ASRF (3, 3>=. 15272 
1120 ASRF (3* 32)=. 12217* ASRF (4*5) =.2018* ASRF (4*7) =. 09727. ASRF <4, 30)=. 0805 
1140 ASRF (4* 32) = . 222. ASRF (5* 6) =1 . 01 45, ASRF (5.7) =. 3354* ASRF (5,32) *. 3954 
1160 ASRF (6* 32) =2 . 59836 , ASRF (7 * 3 0 > =. 1 0 1 23 , ASRF (7 , 32) =. 4939 , ASRF <8, 23) , 

ASRF (9* 32) = . 04729* ASRF (9*24) , ASRF (9, 32) =.0715 
1 1 30 ASRF ( 1 0* 23) * ASRF ( 1 0 * 32) = . 1 0 079* ASRF ( 1 1 , 22) * ASRF ( 1 1 , 32) = . 09727 
1200 ASRF (12,22) .ASRF (12,32)=. 072, ASRF (13*23) .ASRF (13, 32) *. 03857, ASRF <14, 

24) , ASRF < 1 4* 32) =. 05236, ASRF <15, 25) » ASRF < 1 5 , 32) *. 0273 
1220 ASRF (16,26) * ASRF (16* 32)=. 03436* ASRF (17*27) .ASRF (17* 32)=. 06972, 

ASRF (18. 29), ASRF (18,32)=. 06152, ASRF (19,28) .ASTF (19,32)=. 05236, 

ASRF (20,27) * ASRF (20, 32) =. 05903, ASRF <21,26) , ASRF <21, 32)=. 03436, 

ASRF (29, 30) =. 0206, ASRF <30,31)=. 0206, ASRF <30* 32) =. 1296 
1240 LNGTH<1,2)*.0167,LNGTH<2,1)».0375»LNGTH<1*5>=. 0593*LNGTH<5» 1>=. 0917, 
LNGTH <2 , 3) =. 0375 , LN6TH (3,2)=. 0583»LNGTH<2,5>=. 041 7, LNGTH <5, 2) =. 0917, 
LNGTH<3*4)=. 0533.LNGTH <4,3> =. 1 125*LNGTH<3*5> =. 02 08, LNGTH <5, 3)*. 0917, 
LNGTH<4,5>=. 026 , LNGTH <5 , 4) =. 091 7* LNGTH <4 *7) =. 0333, LNGTH <7 ,4) ». 0917, 
LNGTH (5 * 6) =. 09 1 7 , LNGTH (6 * 5) =. 091 7 

1 25 0 LMGTH (3 * 7) = . 2 033 » LMGTH (7 * 5> = . 0333 * LNGTH <4 * 3 0) =. 11 04. LMGTH <3 0 » 4) =. 0375 » 
LMGTH <7*3 0> * . 0333*LN6TH(30*7) =. 0375 
1260 NFM<1)»92.9 
1280 H*15.453»<MFM<1) *♦-. 8) 

1300 HCONV(I,32)=l FDR 1=1 TO 4 
1320 HCONV(I*32)=l FDR 1*8 TO 21 

1 340 HCONV (3*32) * HCOMV (6*32) * HCOMV <7, 32) * HCQNV (30* 32) =2 

1360 HCOMV (8,25) ,HCDMV (9,24) * HCONV <1 0,23) * HCOMV (1 1 *22) , HCOMV < 12, 22> , 

HCONV <13. 23) *HCOMV < 14*24) ,HC0MV<13,25) ,HC0NV(16»26) , HCOMV <17, 27) »H 
1330 HCONV (19, 28) * HCOMV (19*28) *HCOMV(20,27) * HCOMV (21 ,26) , HCOMV <29, 30) , 

HCOMV <30* 31 )=H 

1400 MASFLO (22,31 • ,MASFL0(23.22) ,MASFLO (24,23) .MASFLO (25.24) , 

MASFLO <26 * 25) , MASFLO <27 , 26) .MASFLO <28, 27) .MASFLO <29, 28). MASFLO <31, 29) =1 
1420 I FLUX (3, 1) =45441 , 1FLUX <8, 1) =121 176, IFLUX <20. 1) =15147, 1FLUX<21, 1) *30294, 

1 FLUX (30*1 >*30294 
1440 SSTEMC 1,1) =77, SSTEM (2*1) =652 
1460 HRIM(l) =0 

1480 REAL VOL ( 32) '6.2863,1 0. 234 •• 9. 793,21 .3056*245. 36* 397. 41 * 36.624, 

. 1 36*. 203*. 291*. 251 *.216*. 2557*. 1414*. 079,. 1*. 198*. 1847*. 137, 

. 1676* . 1 05* 1.669*2. 012,1. 197* 1 . 021 * . 7, 1 . 352, 1.216, 1 .216,27.9, 

1,1< 

1500 VOL < I ) »VOL < I ) -d 728 FOR 1=1 TO MM 
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FILE 2 


[ LASTPART ] 


3000 ! CALCULATION OF THE VI EM FACTOR" FOP hXJ SYMMETRIC: NODES 

3020 A ' ni < I • =F‘ T ♦ » R '■ 2*1-1 • + R »: *♦ I 1 * ♦ 3 OR •“ • • R ' c ♦ I > -R • *♦ I - 1 • 1 ♦ ♦* > ♦ 

•: *?<*♦! • -Z*?M-1 • ■ ♦♦*) • FOR I = t TO NH 
3040 FOR 1=1 TO NN 
3 D60 FOR J=I TO NN 

30*0 IF V<I,.h=0 OR V 'I * J • = 4 THEN AQ TO *5*0 
3100 IF V * I • .1 > = 3 THFN 00 TO *5*0 


3120 

IF AES 1 V » I , J > 

. = 

2 THFN AO TO 

3320 


3140 

PH - 1 • =P •: *#.l- 1 

* * 

PH * 2 > =P * *♦ J • . 

PH* 3'=P‘2^ l> 

, RH* 4 * =R * *♦.!- 1* 

3 1 6 0 

PL il ,=p,*#l>. 

RL 

*' 2 ' =P 1 .?♦ I 1 , PL 

• :3 > =P *: 2 ♦ I — 1 ' 

, RL • 4 :> =R * 2 ♦ I - 1 • 

31*0 

S' 1>=Z •*♦.!- 1' 

_ 7 

• *♦1:. • S* 2^=2i 

*♦ J * -7 i *♦ I , 

"■ • 3 ' =2 '* -7 ' *♦ I - 1 > , 


I v4'=7i*# i-l > 

-- 

* *♦ I - 1 * 



3200 

i ll * = 1 + • 'PH* 

LL 

> PL 1 L L i 

* •: S * L L * RL* 

L L ' •'♦♦* > IF RL ML ' :t 0 


FOR L L = t TO 4 

3rc 0 DDF 'LL* * * >■ 'Ll - * * -'"OR *■ ■ i X 'LL' 1 . '2 * ♦ ♦* * - • •: RH »: 1. 1 > Ft ML )'>♦♦*'> • 


IF RL ML) 55 fi FOR LL = 1 TO 4 

3c 40 F . I , ,1 . s . P I A mi • I > - * • * * R I • ♦ ♦£> ♦ ' DDF • 1 • -DDF »:*> < > + 
•: • R - 2* I - 1 ' > ♦ • DDF < 3 :• -DDF <4 :* > > :■ 


3280 IF F «: I * j • 0 THEN F • I , J * = O 

3:300 GO TO *5*0 

3330 IF V •" I , J * = - 2 THEN 11 = J. H = I ELSE 11 = 1* U«J 

3340 I NR= • R. c-# .11 * +R * 2* J1 - 1 ' 2* DTP= * R .;,=■♦ I 1 > +R *2*I 1- p ■ ■ 3 

3 3*0 EiJV'L = . ' 7 • 1 1 -*-7 • : ♦ J1 • - - 7 *?^I 1-1 * + 7 H-1 :• • > 

33*0 XX= J 1 > +7 •. .?♦ 1 1 - 1 ) - •: 7 1 1 * -t-Z • 2* 1 1 — 1 > • > 

3400 K 01 T = • OTR- 1 NR • SAP- * • flTR- 1 NR '■ ♦♦2) + * ' > » R R = □ T P I NR M R = E 0 V L ‘ I N 

342" A=- 1 + • LR ' « B= 1 + • L R P ♦ ♦ 2 • 

344 0 El — * OR 1 * 1 A+2 > ♦ > — ■ 4#PR^R R • • * E 2 = Ai' 0 ' P 1 E • RR ♦h '• 1 i E 3 = £:♦ A _ I NR ■ 1 RR ’* 

34F0 E4=R I#A 2» F5=AC □ : R 1 B A • , E*= 1 • PI#RP> 

34* 0 PH I I=hTNR » * R < £♦ J 1 - 1 • -R • ?♦ .1 1 .» > ■- • 2 -'*♦ 1 1 • -2 .1 1 - 1 • > 

3500 PHI I =ATMP ■ • P 1 * -P '*♦ I 1 - 1 • * • * ?• c‘^I 1 '* -/ ?♦ I 1-1 •• • 

3520 Fill, Jl =LOST^rOSP - PH I I+FHI j - ♦ •: a PR - - > £*,• • F5- < ■ * E1#E2 ' +F3-E4 • 
•*#LP - * ' • > 

3540 FMI1* II' =A SOI. ' 1 1 • ♦F i 1 1 « .11 • A ’ OL • J 1 • IF V < I * I - =-? AND A 3 DL '• J 1 '* 5? 0 
35*0 IF F > I 1 < J 1 * 0 THEN F-Il* J1*=0 
35*0 NEXT J 

3*00 FUM'=F'I*L' FOP L = 1 TO NN 
3*20 FOP \ = 1 TO NN 

3*40 IF V ■ I , y • =4 THEN FiI*F‘=l-iFU#UNV- 
3**0 IF F 1 l M ' 0 THEN F • I , K -=0 

36* 0 F4 i I '=F'IO '♦H" 01.. ■ I • A * OL > \ - IF A' 01. • ► - 5: 0 AND ► I 
3700 NEXT y 
372 0 ne::t i 

374 0 IF VF'RNT = 1 THEN AO TO 37*0 EL “ E AO TO 3 , f*0 0 

37*0 PRINT NON-TERO EL FMFNT I OF THE VIFM FACTOR MATRIX hRF -* 

37*0 FOR 1=1 TO NN 

3*00 FOR .1=1 TO NN 

3*1 0 IF F i I * J > = 0 THFN AO TO ':*4n 

3*2 0 PRINT F ■ : I : , : J: -= r F * 1 * J ■ 

334 0 riF ;T J 
3**0 PRINT 
3**0 ME: T I 

3*00 PRINT r HF("> ON : VMM AT I ON RULE OF F 
3320 PRINT F ♦NN V * 

3 A 40 • COMPLETE THE SFCOND HALF OF THF MATRICFS 

HD J-ASRF-HC ONV 

3^*0 AD J= AD J4-TPM • AD J - * A; RF = A "RF + TRN • A ' RF > * H C 0 N V = H f 0 N V ♦ T R N - nr ONV • 


3**0 • 

4000 • CONVERT INITIhI TFMP IN DEG. F TO DF A . AF f DLU TF 

4 1 1 2 i’i TFM=TEM+ 4*0#MNV. 

4 04 0 • 

4 0*0 * COMPILE F CIRCUMFLEX MATRICES FOR I NR RAPED RADIATION AMD 10LAP 
40*0 IF ‘ h I PR- 0 THFN AO TO 42*0 
4 1 00 RH0= MIHV-EIP ' • PHD - T= • UNV-ESOL > 

4120 FI I , J> =V D ' I « J*- ' F ■ T , J* ♦RHO* I* * FOR T = 1 TO NN FOR. J=1 TO NLJ 
4140 F 1 " • I , J * =► D • 1 « J ■ - * F 1 I « J • ♦RHn r T •' I .* ) FOP T = 1TD NN FOR ' J= 1 T0 - NN 
4 1 * 0 F CF=F • I NV • F 1 « , FC F I T = F ♦ I NV 'FI S ' ‘ . ! * ** • ' ■ ? 

41*0! DF TERM INF INFRARED EXCHANGE COFFFIjSIEHT 

4*00 IREXC . I, j-. = ’ IGMh^R IP. I > ♦EIR- J* ♦ASOL. • 1 ■ ♦FCF * I- J> FOR 1 = 1 TO NN 

FOR J=1 TO NN 


42*0 IRFXC 1 = IRFXf - • * I RE: :C ♦MNM > MUL ► D ' 


ORIGINAL ? A6E J? 

OF *00R 
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4340 ! 

4360 IF CONDUC-O THEN GO TO 4360 

4380 ! COMPUTE CONDUCTION COEFFICIENTS FOP ADJACENT SOLID NODES 
4300 CnNDC < I . J > *ASPF < I • .h ♦KCOND < I ) 4KCOND < . » x < CKCOND < J> *L N6TH < 1 » J > > ♦ 
<(cCONDa>4LNGTHCJ» I>)» IF ADja»J>»l FOP I»1 TO NN 

FOP J*1 TO NN 


4330 CONDC“CONDC- < <CONDC*UNM> MUL KD) 


4340 I 

4360 K*1.SSK*1.ENEX*0.SFIN«0. INTEP»0.EFF«0 

4380 CAPV*DEN MUL VOL MUL SPHT 

4400 IF S?T*t THEN CLOCK*HPINOO 

4430 IF IONDIN THEN GO TO 4460 

4440 IF CL0CK>«HPINOO THEN 60SUD 5740 

4460 IF CLOCK-HPOUT THEN PPINT s'TINE <HP> ' « TAB < 13> J 'AQC ENGY EXTP<BTU> ' * 
TAB<34>t'ACC SOL ENGY IN<BTU> 'iTABC58> * 'PECVP EEF'\HPQUT»TAB<19> « 
FNEX* T8B <'40') ! SEINt TAB ^GO’) tEFF 

4430 IF CLOCF-HPOUT THFN PPINT \'NDDE NO' * TAB < 13> J 'NODE TEMP.F'x 
k.'K ! TAB <13>: TEM 000-460 FOP Kkr*l TO NN 
4500 IF CLOCk=HPOUT THEN HPOUT»HPOUT+DTOUT 
4530 IF HR0UT<CL0CK+DT8U AND HROUT>CLOCK 

THEN TEM1»TEM. CLOCK 1 -CLOCK . INTEP-I »ENEXI*ENEX. SEIN1-SEIN 
4540 ! COMPUTE CHANGE IN TEMPERATURE AT EACH NODE DUPING THE 

TIME INTERVAL DTAIJ AND CALCULATE NEW TEMP 
4560 IF SST=1 AND SSK<*NDIN THEN GO TO 5060 ELSF IF SST*i AND 
SSIONDIN THEN GO TO 5300 
4580 ! TRANSIENT SOLUTION 

4600 IF 2CNI=0 THEN GO TO 4640 

4630 ?CN ■" I > * 1 IF PTAUVaxMINDTU AND SOPS IN < I > *0 FOP I“1 TO NN 
4640 PES«<IPEXC1»<TEM MUl TEM MUL TEM MUl TEM) > +POWP+SOLCON+ 
i (Ft n+CONDC+CONVO *TEM> 

4660 TEMNFW=TEM+<DTAU»PES div capv> 

4680 GO TO 4930 IF 2CNl=0 

4700 DTEM=TEM 

4730 FOP JJ*1 TO MNIT 

4740 RES* 1 1 PEXC 1 ♦ <DTEM MUL DTEM MUL DTEM MUL DTEM) > +POWP+SOLCON4 
i (FLO+CON DC +CONVC ) *DTEM) 

4760 JAC Cl* J>*<4*IPEXC1 >ri , J) ♦ «'DTFM < J) ♦♦3’* > ♦FLO < I • J> +C0NDC <1 « J> ♦ 


CONVC I • J > FOP 1*1 TO NN FOP J*1 TO NN 
4780 FOP 1*1 TO NN 
48 0 0 RES < I > * ft IF ?CN < I > * 0 

4330 JAC < I » L > “0 IF ?CN<1>»0 OP ?CN<i>«0 FOP L*1 TO NN 
4340 JAC < l » I ) * 1 IF 7CN ( I > *0 
4360 NEXT I 

4880 DTEMsDTEM- • I MV* JAC) ♦PES) 

4900 NEXT JJ 

4930 FOP 1 = 1 TO NN 

4940 TEM < I > =DTFM • I » IF 7CN ‘ I > *1 

4960 TEM < I ' *TFMNFW ( I > IF S0PSIN*I>*0 AND ZCN<I)»0 
4980 NEXT I 
5 0 0 0 C L DC: ¥ *C t DC ¥ ♦ p T HU 
5020 GO TO 5?60 

5040 S STEPpV \THTF ^01 LIT I ON 

5060 FOP J.l*i TD MNJT 


5030 


5100 


PEWIRFXCUlTEN MUl. TFM MUl TEM MUl TFM) > ♦POMP ♦m CON* < <FLO* 
CONPC+CDNVO *TEM> 

JflC •" 1 1 .1 > * <4# I PEXf: l *: l • . 1 ) ♦ < TFM t J ) :< > ♦FLO < JO* ♦CCINOC • I • J * ♦ 


CONVC <1% .¥> FOP 1*1 TO NN FOP »*l TD NN 
5120 FOP I«l JO NN 
5140 PEC • P * 0 IF -nP' lNc I>»l 


5160 I *L>*0 IF SOPSIN •: I > *1 OP SOPSIN*t>*l FOP 1*1 TO NN 


5160 J*C< I • J > ** t IF *OP!lN<J>*i 
5200 NEXT J. ' •/ 

5330 TEM-TEMr^I^*: lftC>'*PES» 

534ft NEXT JJ 

5360 PES* < I PEXC 1 ♦ ''TEM MIJL *TEM MIJL TEM Mill TEM) ) ♦PtlMP+SOLCON-) <' <Fl Q+ 
CONDC+CONVC* ♦TEMC 

5330 PRINT v' NODE NO 1 PAB ‘ I3> « NODE RESIDUE • BTUH' s FKl TAB < I3» J PES <KK> 
FDP KI *1 TO NN 

53ft0 CLDC* *CI OCK+riTAU. SSF-SSK + 1 
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5320 i 

5340 1 CALCUL AT ION DF ACCUMULATED ENERGY 

5360 IF SST-0 AND FNDTN«0 ThFN ENEX-ENEX* <DTAU4MFM <K-1* ♦ SPHT <FNDIN> ♦ 
< TFM CFNDOT > -TFM <FNDTN> > > 

5330 IF SST-1 THEN ENEX-ENEX+RES <FNDIN> 

5400 SEIN-SEIN-kDTAU-CUNVrCFL.UX MUL ASOL>V* 

5420 EFF-ENEX'SEIN IF SEIN»0 
5440 IF INTER-1 THEN GOSUB 5520 
5460 IF CLOCK <=FTIME THEN GO TO 4420 
5430 END 


5500 HINEAP INTERPOLATION SUBROUTINE FOR PRINT OUT 
5520 TEMIN-TEM1 ♦ cHROUT-CL OCK 1 > ✓ ^CLOCK-CLOCK 1 > ♦ CTEM-TFMl* 

554 0 ENEX l N-FNFX 1 ♦ < f HROUT-CL OCK 1 > ✓ CCLOCK-CI OCX 1 > > ♦ < ENEX-ENEX 1 > 

556 0 SEININ-SFIN1+ < CHROUT -CL OCK 1 > CCLOCK-CL OCK1 > > ♦ <"SEIN-SEIN1 * 

5530 FFF-ENEXIN'SEINTN IF SEININ«0 

5600 PRINT \'TIME<HR*'»TAB<12>«'ACC ENGY EXTR CBTU> ' I TAB <34> t 

'ACC SOL ENGY INCBTU* '«TAB<58> J 'RECVR EFF '^HROUTt TAB* 19* : 
ENEXIN: TAB *40* s SEININS TAB <60* JEFF 

5620 PRINT v'NDDF NO' « TAB <12* : NODE TEMP.F' vKKi TAB < 12> « TEMIN • KK > -460 
FOR KK-l TO NN 

5640 HROUT -HROUT +DTOUT • INTER-0 
5660 RETURN 
5630 • 

5700 ! 

5720 1 TIMF DEPFNDENT VARIABLES AND STABILITY SUBROUTINE 
5740 TFM <AMBNN> -S8TFM •: t « K > +46 0 « TEM •: FND I N> *SSTFM*2«K >*460 
5760 7CNI-0 

5780 ZCN* I> =0 FOR 1=1 TO NN 
5800 IF SAIRR-0 THFN GO TO 5930 
5820 FI UX< J>»IFLUX< J.K* FOR J-l TO NN 
5340 ! 

5360 ! COMPUTF CONSTANT TERM OF FINAl VFCTOR EGUATION 

5380 SOL CON-ASOl Mill <<ESDL MUL • FCFST4 LRHTIST MUL FLUX'* >*♦ 

• FSOL MUL FLUX* > 

5900 L IRFXC <1 » I* -IREXC < I . J* ♦ < *TEM< J» ♦♦2* ♦ <fTFM<I > *42* >»<TEM< .1* ♦TEM 1 1 > > 
FOR 1-1 TO NN FOR .1-1 Tn NN 
5920 l IREXC-l IREXC-< a IREXC^'INM* MUL KB* 

5940 1 

5960 ! DETERMINE CONVECTION AND FLOW MATRICES 
5980 IF CONVFC-O THEN GO TO 6040 

6000 CONVCa. J*-HCONVa« J»RASRFCI, J> IF fllll('I« l*»2 FOR 1 = 1 TO NN 

FOR J=1 TO NN 

6020 CDNVC-CONVC- < <'CONVC«UNN> MUL KD» 

6040 IF FLOW-O THFN GO TO 6100 

6.060 FLO <1 » J * -MASFLO < I * J * ♦SPHT < I ■ -MFM < k > FOR 1*1 TO NN FOR .1=1 TO NN 
6030 FI 0=FI_0-< CFI n-IINM* Mill KB* 

6100 GO TO 6280 IF SST-1 
6120 ! 

6140 ! COMPUTE MIN TIMF INTERVAL DTAH FOR STABLE TRANSIENT SOLUTION 
6160 OT AUV < I * — 1 ♦CAPV < I > ' <FI. D < I . I > ♦CONDC • I . I * ♦C'ONVC < I • I > H I PFXC < I . I ■ • 
FOR 1-1 TO NN 

6180 PRINT V NODE NO' « TAB <12* t 'NODAL MIN DT.HP '\KKl TAB ■ 12 » i DTAUV • ► K * 
FOR KK-t TO NN 
6200 DTAU-1. 

6220 DTAU-DTAIJVCI L * IF DTAUV <1L > <DTAU AND SORBIN <LL > -0 FOP LL = 1 TO NN 

6240 IF DTAU MINDTU THFN DTAIJ-M INDTU. ZCNI-1 

6260 PRINT MIN TIME INCREMENT IN SECONDS IS • J • DTAU*3600> 

6230 IF SST-1 AND K>1<-NBIN THEN DTAU- ‘HP I N *K ♦ 1 > -MR I N <:K * > ELSE IF 
SST«1 THFN DTAU-1 
6300 K-K>1 
6320 RETURN 
> 
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